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Abstract 


This report describes the algorithms used in the post-Block II fringe-fitting software 
called “Fit.” The steps needed to derive the VLBI charged-particle corrected group delay, 
phase delay rate, and phase delay (the latter without resolving cycle ambiguities) are 
presented beginning with the set of complex fringe phasors as a function of observation 
frequency and time. The set of complex phasors is obtained from the JPL/CIT Block II 
correlator. The output of Fit is the set of charged-particle corrected observables (along 
with ancillary information) in a form amenable to the software program “Modest.” 
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Section 1 


Introduction 

This report presents the theoretical basis for the post-Block II fringe fitting data anal- 
ysis used in very long baseline interferometry (VLBI). The overall goal of this analysis is to 
extract unbiased physical observables for each observation scan, namely the charged par- 
ticle corrected bandwidth synthesis (BWS) delay, phase delay rate and phase delay (when 
cycle ambiguity resolution is possible), directly from the Block II correlation coefficients. 
These observables are typically combined with those from other scans and passed to the 
fitting program Modest' 1 ! (formerly called Masterfit) for further analysis. 

By design, much of the material in this report is independent of a specific software 
implementation. However, where choices in computer memory allocation or processing 
time influence the analysis, this report will refer specifically to the program Fit and its 
solutions to the various data analysis problems. As experience with Mark III data increases 
and as better computational resources become available, detailed improvements will be 
developed and implemented in Fit. Although some details of the implementation may 
change as future improvements are made, the majority of the material presented here 
should remain current. 

The Block II processor is a VLBI correlator developed at JPL and the California 
Institute of Technology. A detailed description of this processor can be found in references 
2 and 3. The Block II is capable of correlating both Mark II and Mark III data formats! 3 ', 
although this report will concentrate on the Mark III data type because its higher data 
rate makes it much more desirable for high accuracy work (Fit is able to process Mark II 
data and much of the material presented in this report is valid for both data types). 

Some aspects of post-Block II data analysis have been described in reference 2; however, 
at the time of its writing, the Block II was expected to include a hardware device named 
Tensor. This device was designed to perform the Fast Fourier Transforms needed in the 
data analysis described in sections 6 through 11 in this report and section 13 in the 
reference. Work on Tensor was abandoned before its completion, leaving the post-Tensor 
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software, named Restore, without an input route, except for a software emulator of the 
Tensor hardware. This software emulator was designed for debugging Tensor during its 
construction, and not for data analysis. Thus, except for the software emulator of the 
Tensor hardware in combination with Restore, no route from the Block II processor to 
VLBI observables existed. 

Because Tensor’s function would eventually have to be performed by post-correlation 
software, the opportunity to drop many restrictions and approximations required for such 
a hardware device was exploited. This report presents these new analysis steps, along 
with all those needed to obtain the final observables from the Block II output, and their 
theoretical justification. The material presented here can be roughly divided into three 
tasks: a) obtaining an a priori estimate of residual delay, phase delay rate and phase 
for each observation band, b) calculating accurate corrections to the a priori estimates 
using phase tracking and linear least-squares parameter estimation, and c) combining these 
estimates with the correlator model and applying charged-particle corrections to obtain 
the final observables. Each of these tasks is divided into several topics which are covered 
in detail by each section of this report. Sections 2-5 cover some preliminary calculations 
necessary to the analysis, sections 6-11 discuss the a priori calculation, sections 12-15 
are concerned with parameter estimation, sections 16-18 pertain to the final observable 
estimates and section 19 covers the Modest interface. Appendices A and B describe specific 
calculations in more detail than that found in the text, appendix C gives the mathematical 
definitions of special variables used extensively in the least-squares analysis and charged- 
particle calibration, and appendix D presents a complete list of variable definitions. 

The Phase Sign Convention Used in this Report 

Because observables extracted from the Block II output are passed to Modest, the sign 
convention of both must be treated explicitly. The Block II’s phase sign convention is sta- 
tion 1 phase minus station 2 phase, while Modest’s is the reverse. Most expressions in this 
report are valid in either sign convention as long as all quantities are in the same system. 
However, the Block II convention is implicitly assumed for all residual observables because 
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these are derived more directly from the correlator output^ total observables, because they 
are passed to Modest, are assumed to be in the Modest system. This separation of the 
two sign conventions corresponds to that used by Fit. In this report, the Block II system 
is used up through section 15 and the Modest convention is used in sections 16 to 20. 
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Section 2 


Correlation Coefficients and Error Analysis 

This section presents expressions for the complex correlation coefficients and their 
errors as a function of the number of Block II accumulation counts: the Block II counts 
are the raw inputs to Fit. Expressions relating correlation coefficient errors to phase and 
amplitude errors are also given. Finally, the variance of instrumental phase calibration 
phasors and their effect on the calibrated cross-correlation data are derived. 

The Frequency Domain Correlation Coefficient and its Variance 

The delayed bitstreams from each station at the correlator are multiplied together and 
counterrotated! 2 ! j n quadrature for a number of lags, typically eight. A Fourier Transform 
(FT) is performed to transform these lag-domain correlation coefficients to the frequency 
domain. These final frequency domain correlation coefficients are the primary correlator 
output passed to Fit. The frequency domain correlation coefficients for bin 6, can 

be written 

i n l~ 1 ■ bl 

m = 1211 

L 1=0 

where the sum over lag shows the FT explicitly, pi(t) is the complex correlation coefficient 
for lag l after fringe counterrotation at time t and Nl is the number of lags. Both b and / 
are integers ranging from 0 to Nl — 1. For the typical case of eight lags and a sample rate 
of 4 MHz, each of the eight bins has a nominal width of 4 MHz /Nl = 0.5 MHz. Because 
the channels are single-sideband, four of the eight bins contain noise and are ignored by 
Fit. After expanding the real and imaginary parts of [2.1], the variance of pb(t) can be 
written 

Var[Re />,(()] = — p-j £ [var[Re ?,(!)] cos 2 (^) + Var[Im ?,(()] sin 2 (^)] [2.2] 

L 1=0 L L 

N l —\ 

Var[Imp 6 (i)] = ^2 J2 [Var[Im £,(*)] cos 2 (^) + Var[Repj(i)] sin 2 (^)j . 
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This step assumes each real and imaginary component is statistically independent of every 
other component, real or imaginary, as derived in Reference 4 (this has been tested by 
measuring the autocorrelation of the Block II bitstreams! 5 !). The lag domain correlation 
coefficient after counterrotation (also called lobe rotation), pi(t), can be expressed in terms 
of the correlation coefficient before counterrotation as 

V 

E i * i J fTO 

[2.3] 

3 = 1 

where 0i{t b ) is lag Vs single-bit correlation coefficient before lobe rotation, <f>t m is the model 
fringe phase, t b is the time tag of the jth bit in the sum and V is the total number of bits 
in the sum, typically around 8 x 10 6 for two-second integrations. Assuming the observed 
bandwidth is not oversampled and no DC bias exists (so that each bit is statistically 
independent and each component of the correlation coefficients has zero mean), the variance 
of the real and imaginary parts of pi(t) can be written 

v 

Var[R ep,(t)] = Var [&(*$)] cos 2 [2.4] 

y= i 

V 

Var[Imp,(t)] = ^ Var^tjjjsin 2 <^ m . 

;=i 


Equations [2.4] can be evaluated from the variance of the like and unlike-sign counts 
which are binomially distributed. For the digital three-level lobe rotator used in the Block 
II correlator, the cos 2 <f>t and sin 2 <f)^ m terms will be +1 when not blanked and zero when 
blanked (see references 2 and 6 for more detail on lobe rotation). For an interval of V bits, 
the real and imaginary components of equation [2.3] become 


Re pi(t) = ^ ~ = ( 2 ^x + - V) 

im Pi(t) = ^ ~ 


[2.5] 
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where V = L x + U x + B x = L y + U y + B y and L x[y) and U x[y) are the number of like and 
unlike-sign counts of the real (imaginary) component of p t (commonly referred to as the 
cosine (sine) accumulation) and B x{y) is the number of blanked bits in the mterval. From 
[2.5], the variance of the x(y) component of pi is 4 /V 2 times the variance of £*(„)• Using 
the binomial distribution’s expression for the unbiased variance estimate of the number of 
successes in V - B x trials, given the actual number of successes is L z , gives 


Var[Repi(t)] 


4 

V 2 


( 


V - B x 

V - B x - 1 




Lx 

V - B x 




[ 2 . 6 ] 


with a similar expression for the variance of the imaginary component of p,. Solving the 
first equation in [2.5] for L x and substituting leads to 


Var[Re pj(t)| = 


V -B x 

V 2 


(S) Mi^)Va(‘>i’]- 


(2.7] 


Using V — 5j( y ) > 1. pi < 1 and the imaginary component analog to [2.7] gives, to very 
good approximation 

V — B 

Var[Re p/(f)] ^ Var[Imp/(t)] ^ , [ 2 - 8 ] 

where B is the average of the blanking counts for the sine and cosine accumulations; for 
typical 2-second integration intervals and several KHz fringe rates, the two blanking counts 
are nearly identical. Finally, substituting this into equations [2.2] leads to 

a 2 = Var[Rep 6 (*)] ^ Var[Imp 6 (t)] c; , [2.9] 


where a 2 is the variance of either component of p b . This shows the system noise associated 
with a given fringe phasor depends only on the number of processed bits and number of 
lags (given the above approximations). For the Block II lobe rotator, B « U/4 so 


a 


2 


3 

“ 4N l V' 


[2.10] 
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Expressing the Correlation Coefficients and their Variances in terms of Corre- 
lator Output 


For each dump interval, the Block II correlator passes the user the like-sign and blank- 
ing counts for the sine and cosine accumulations along with the invalid counts. The invalid 
counts were ignored in the above analysis but can be included by interpreting V as the 
number of valid counts: the total number of bits, N, minus the invalid count, 7. The real 
component of the correlation coefficient can be written as (see [2.5]) 


Rep/ 


L x - U x 

N -I 

2Lx — TV + 7 + J? z 
N-I 


[ 2 . 11 ] 


with a similar expression for the imaginary component. The Block II correlator hardware 
has a flip-flop circuit before each accumulator which results in the counts sent to the user 
being only half the expected value. Therefore, L x{y) = L x{y) f 2, B x{y) = B x{y) /2 and 
I I where the hat ( ) represents quantities passed to the user by the correlator (the 
correlator accumulator flip-flop circuits alternately begin with 0 or 1 to avoid inducing a 
bias, on average). Thus, 


Re pi = ^ Af/4 + 1/2 + B x / 2 

Nf 4 - 7/2 

L x — Z x 

~ M ’ 


[ 2 . 12 ] 


where = N/4 - B x / 2 - 1/2 is a DC bias and M = N/4 — 1/2 is a normalization 
factor ( Z x is a DC bias in the sense that the Block II hardware accumulators can only 
increment, and must therefore be offset by Z x , so the cosine accumulator contains L x — Z x , 
for example). The expression for the imaginary component of p, is similar. After the lag 
FT, given by [2.1], the fringe phasor is 


Pb 


( Pz,Py ^ (-N/4 — //2* N/4 — 7/2^ (m’m)’ 


[2.13] 
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where S and C are the like-sign sine and cosine counts from the correlator (the Block II 
can process data in either “lag” or “bin” mode; the sine and cosine counts passed to the 
user are given by L x and L y in equation [2.12] when in lag mode, and C and S in [2.13] 
when in bin mode, which is the mode assumed in this report; in both cases, C and S 
are integers, because, when in bin mode, [2.1] is calculated using integer arithmetic). The 
error on each component, obtained from [2.9] above, is 

2 N-I-B V-2 B M-B/2 f214 | 

a - Nl ( N - I ) 2 " N l V 2 4 N L M 2 

where V = N — 2/ = 4Af. The software program Fit uses [2.13] to calculate the correlation 
coefficients and [2.14] to calculate the distribution variance, with B equal to the numerical 
average of the sine and cosine blanking accumulations, as noted above. 

Because of the large number of bits processed by the Block II for each phasor in a typical 
two-second dump, usually about 8 x 10 6 , and the inherently small correlation coefficient, 
the central limit theorem implies the probability distribution of either component of p b is 
normal to very high accuracy. The real and imaginary components of p b (t) are therefore 
independent, as proved in Reference 4, and have a Gaussian probability distribution. The 
means of the real and imaginary components, p x and p y , are the components of p b given 
in equation [2.13] and cr 2 is the variance of either component as defined in [2.14]. 

Phase and Amplitude Distributions and Variances 

Given the Gaussian probability distributions for the real and imaginary part of the 
fringe phasor, with means and variances given by [2.13] and [2.14] above, it is possible 
to change variables to derive the distribution for the phase and amplitude. Appendix A 
contains a complete derivation of the phase and amplitude probability distributions which 
are summarized below. Defining <p = arctan(p y ,p z ) and A 2 = p\ + p 2 for a given bin, the 
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phase distribution, mean phase, {</>}, and phase variance, are given by 


1 -ft 2 

P(4) = 7T e 2,72 

Z 7T 

W = ^ 


I*R ,, aR\ SL co*‘[*-+ a ) , R cos(<j> - <t> R ) 

l + — cos{4> — 4> R )e 2 <r 2 M + erf - ' 


/ i_lR i 


<4>-<t> R ) \] 

V2a )\ 


o 7T 2 nfti 
°<t> ~ ~^~ e 2a + 


R 


y/2fta 


/' 


-R 2 »in 2 9 


d0 6 i cos 0 e 2 ^ ( 1 _l_ er f 


R cos 0 s 
\/2 a ‘ 


[2.15a] 

[2.156] 

[2.15c] 


where 4> R and R are the polar coordinates of the parent probability distribution’s center 
(see below). The amplitude distribution, mean, and variance are given by 



where I 0 and h are modified Bessel functions of zero and first order. Table 2.1 shows some 
of these expressions for the no signal, weak signal and strong signal limits. Figure 2.1 
shows a# plotted as a function of R/a along with the weak and strong signal limits. Note 
that R/a is defined to be the signal to noise ratio (SNR) for a given phasor. 

In the above expressions, a careful distinction between R and A should be made. 
The parent probability distribution for p b is a two-dimensional Gaussian centered at polar 
coordinates (. R,4 > R ). In the no-signal case, R is zero. The amplitude of a given phasor, A, 
drawn from the parent distribution will always be non-negative, even in the no-signal case. 
Thus A is a biased estimate of R as seen in [2.16b] or the first row of Table 2.1. Figure 2.2 
shows (A) /a as a function of R/a. 
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| EXPRESSIONS FOR (A), a a AND o $ IN LIMITING CASES 

Statistic 

White Noise 

Weak Signal 

Strong Signal 

(A) 


a VT 

+ &) 

R U + ^ 2 ) 

<*A 


< V 2- 2 

°\j^‘ — f (i + ^ 2 ) 




7T 

v/3 


a/R 


Table 2.1 Expressions for the mean amplitude, amplitude sigma and phase 
sigma in the no signal, weak signal and strong signal cases. R 
is defined in the text and o is given by [ 2 . 14 ]. 


Error Analysis of Instrumental Phase Calibration Data 

Instrumental phase calibration tones are injected near the front of the station elec- 
tronics. These tones traverse the same instrumentation as the actual data so phase or 
delay instabilities induced by the electronics at each station can be analyzed (see section 
4 of this report, and references 2, 5 and 7 for more detail). The Block II lobe rotates the 
bitstream from each station in quadrature by a model of the tone phase, for each tone. The 
counterrotated bitstreams are summed over many bits, typically 8 x 10 6 , to raise the SNR. 
This process is similar to the cross-correlation counterrotation except instead of a 3-level 
digital lobe rotator, the tone lobe rotator uses 128 levels. Following similar arguments as 
above (equations [2.3] through [2.10]), the variance of station t’s tone phasor is 

o\ = Var[Rept f (t)] ^ Var[Imp ti (t)] st — [2.17] 

where p ti is station t’s tone phasor and V t . is the number of valid tone counts. This result 
is analogous to [2.10] above; N L is 1 because there is only one “lag” for tone processing 
and the £ compared to \ is the result of better modeling of the sine and cosine functions 
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Fig. 2.1 Plot of a^ as a function of R/a (solid). The strong signal limit 
(dashed) and weak signal limit (dotted) are also shown. 



Fig. 2.2 Plot of (A) /a as a function of R/a (solid). The strong signal 
limit (dashed), weak signal limit (dotted) and (A) = R (dot- 
dashed) lines are also shown. 






for tone lobe rotation (these factors come from averaging sin and cos 2 <^ m for the lobe 
rotation model in the cross-correlation data and similar factors for the tone models). 

When calculating [2.17] from the correlator output, multiplexing of the tone lobe ro- 
tator must be considered. The tone lobe rotator is time multiplexed over 4 tones (no 
multiplexing is done for the cross-correlation lags), resulting in hardware counts from the 
correlator being a factor of 4 smaller than otherwise expected. The analogous expression 
to [2.14] for each station is 


where a ? is the variance of either tone phasor component, /«,. is the station tone invalid 
count and the 8 results from a factor of 2 for the flip-flop circuit and 4 from multiplexing. 


One option in Fit is to have each cross-correlation phasor counterrotated by the ap- 
propriate tone phasor from either one or both stations (section 4 covers this procedure 
in more detail). After calibrating the cross-correlation phasor by counterrotation, the 
resulting phasor has variance given by (see [4.3]) 


o 


2 


a 2 + 



2 


< + 




[2.19] 


where the 1 and 2 refer to the two stations, A is the cross-correlation amplitude, A ti is the 
station i tone amplitude, and <7 2 and o ^ are given by [2.14] and [2.18]. If the tone data 
are used for calibration, the replacement indicated by equation [2.19] should be assumed 
implicitly in place of [2.14], 
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Section 3 


Band Reference Time and Frequency 


A reference frequency for each observation band and an overall scan reference time are 
necessary at many stages of the analysis. 'The expressions for these references are sum- 
marized here because they are needed immediately for the a priori calculation; derivations 
are postponed to later sections where necessary supporting material is given (see sections 
14 and 17). Additional reference frequencies are required for each band and each observ- 
able (BWS delay, phase delay rate and phase delay) for the charged particle calibration 
described in section 17; expressions and derivations for these frequencies can be found 

there. 

The scan reference time, to, is the only time reference necessary in the analysis and is 
given by (see [14.11] and Appendix C) 


to 


(*) 

E 


Ujtj 

t; 


w 

a\ 

0 


Xs 

Xi’ 


[3.1[ 


where t runs over all frequency bins in the highest frequency band (assumed to be X band), 
j mns over all time points, f t , is bin t’s observation frequency, is the phase variance 
given by equation [2.15c], and X 5 and X\ are the X-band counterparts to the expressions 
given by [C.l[. Missing or bad data points are implicitly excluded from the sums. 

The X-band reference frequency is given by (see [14.12] and Appendix C) 

X19 — toXie _ Xtz [o 0 i 

X09 — 2 * 0^06 + *0^02 

with a similar expression for other observation bands (u>os ^ the S-band reference frequency 
and u>o is a generic band’s reference frequency, for example), and the undefined quantities 
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are given by [C.3] and [C.10]. If t 0 is defined as in [3.1], this expression for the highest 
frequency band simplifies to 


U 0X 


(*) 

E uj 

al 

tj 9 ij 

W . 

tj 9t 3 


Xi 

X Q ' 


[3.3] 


The above expressions are evaluated early in the analysis to give a preliminary t 0 and 
w o used in the a priori calculation. After the phasor time integration, described in section 
12, these expressions are reevaluated because some phasors may have been combined or 
eliminated, changing the summations. The new t 0 and w 0 ’s calculated there are the final 
values for these references. The change in these references caused by their reevaluation 
induces implicit shifts in the a priori delay, phase delay rate and phase which must be 
accounted for, as described in the last subsection of section 12. 
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Section 4 


Instrumental Phase Calibration 

The instrumental phase calibration tones injected near the front-end of the station 
electronics are processed by the Block II correlator by counterrotating the bitstreams from 
each station by a tone phase modelM. Residual phases about this model are passed to the 
user as like-sign sine and cosine sums. The counterrotated tone phasors for station t are 
calculated similarly to [2.12] since there is no lag FT: 

_ _ 64 f Ctj -M ti Stj ~ M ti \ [ 4> j] 

pt ‘ ” 63 V M ti ’ M ti ) ’ 

where Mt = iV/16 — 7^/2 is the tone normalization, Ct { , and 7^ are the sine, cosine and 
invalid tone counts directly from the Block II, and N is the total number of bits processed. 
The factor of 64 /63 is the result of hardware accumulator offsets J 5 1 The variance of each 
phasor component is given by [2.18]. 

The cross-correlation phases are corrected by counterrotating by the phase of the coun- 
terrotated tone phasors, for each bin and integration interval. The Block II is capable of 
performing this correction automatically every 25 milliseconds. The disadvantage of this 
operation is that the user cannot reconstruct the counterrotation on that time scale for 
post-correlation analysis; the Block II only passes the tone counts each dump interval, 
typically every two seconds. To have more complete control over the analysis, Fit allows 
the user to apply the two-second tone phasors to the cross-correlation data in software 
rather than during correlation. Actually, Fit allows the processing of any combination of 
the three available data types for each baseline: the cross-correlation data and the tone 
data from each of the two stations. This also allows detailed study of the tone data alone 
before its application to the cross-correlation data. The equations expressing the combined 
phasors’ components and variances in terms of the cross-correlation and tone data are 


4> *— <t> — <t>ti + 4>t 2 

[4.2a] 

A ■*— A 

[4.261 

o + — A v 

(42cl 
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where 


l / 2 = 



[ 4 . 3 ] 


{A,<f>) and ( A t it <f>ti) are the polar coordinates of the cross-correlation and station t tone 
phasors, and a 2 and o], are the cross-correlation and tone variances. If the tone data 
are used for calibration, the replacements indicated by equations [4.2] should be assumed 
implicitly in place of the uncalibrated variables. Equations [4.2a], [4.2c] and [4.3] are 
applicable for any combination of cross-correlation and tone data if one sets <f> = a — 0 
when cross correlation data is not desired in the combination, and (f>t i = Of = 0 when 
station i tone data is not wanted. Equation [4.2b] is applicable only when cross-correlation 
data are in the combination; otherwise A <- y/A^A^ when tone data from both stations 
is selected or A <- A u when just station i is used. These calibrated phasors are now 
processed without regard to the data type(s) in the combination except for the phasor 
time integration, explained in section 12, which involves summing the cross-correlation 
and station data before combining them as described here. 


4-2 



Section 5 


Manual (Self) Phase Calibration 

Each of the Mark III video converters introduces an arbitrary phase in the channel it 
records. In addition, instrumental phase non-linearities across the 2 MHz bandpass induce 
phase offsets for each bin in a given channel. These phase offsets can be calibrated with 
instrumental phase calibration discussed in the previous section or they can be calibrated 
by the method presented here, or both methods may be used. If these phase offsets were 
constant in time and data drop-out effects were negligible, there would be no need for this 
calibration; it will be shown below that drop-out effects are large, making this calibration 

almost always required. 

Figure 5.1 shows a sample time-averaged phase (shown with error bars) as a function of 
observation frequency. For simplicity, phase non-linearities across the 2 MHz bandpass are 
not included here. Three S-band channels with four bins each are shown in this simplified 
example. Because the phase is determined only modulo a cycle, several possible cycle 
assignments are also shown. The residual (about the correlator model) BWS delay is the 
slope of the line connecting these phase points. The dashed line shows the best fit to these 
data, without channel phase corrections, giving a residual BWS delay that is different than 
the slope obtained from the four bins in any single channel, shown by the dotted lines (the 
slope of the dotted lines is referred to as the bit-stream alignment (BSA) delay for the 
given channel). The large bias between the BWS and average BSA delays is caused by 
the channel phase offsets, which can be seen in Figure 5.1 as vertical offsets between the 

dotted lines. 

The residual phases about the dashed line in Figure 5.1 show systematic upward trends 
for the four bins in each channel. If these trends are constant in time, the residual delays 
for all scans will have a constant, clock-like term added to them. In practice, the trends are 
not constant in time due to data drop-outs, resulting in significant delay estimate errors. 
This can be seen dramatically by assuming the data in the first channel were lost. The 
best fit line in this case is shown by the dot-dashed line, and results in a 250 nanosecond 



Fig. 5.1 Plot of phase, shown with error bars, as a function of bin fre- 
quency for several cycles. In this simplified example, the dashed 
line shows the best fit line (BWS delay) for all three channels, 
while the dot-dashed line shows the best fit line for the two high- 
est frequency channels. The dotted lines show the BSA delays 
for the three channels; the differences in vertical offset between 
these lines correspond to the channel phase offsets. 


change in residual delay from the slope of the dashed line. Even the small fraction of 
data drop-outs caused by tape playback at the correlator can cause changes in the phase 
variances which translate into significant biases in residual delay; a drop-out error causing 
the fit BWS delay slope across 100 MHz to shift one degree at one end of the band leads to 
a 28 picosecond error. The sensitivity to data drop-outs can be minimized if the residual 
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phases are all as close to zero as possible. 

The program Fit minimizes the phase residuals by solving for a constant phase offset 
for each bin in a specified scan. This removes the bandpass phase non-linearity in the bins 
of a single channel as well as the channel phase offsets. A straight line with arbitrary slope 
can be drawn on Figure 5.1 and the phase differences between the data points and this 
line can be taken as the set of bin phase offsets; subtracting these offsets from the data 
will cause the phase points to lie on the line. Because the slope of the line is arbitrary, 
this method introduces an arbitrary clock-like term in the final delay estimate for each 
observation band. By using the specified scan to determine the bin phase offsets, and 
using these offsets for all scans in the experiment, the arbitrary clock for each band is 
equal for all scans. If a high SNR scan is used to obtain the bin phase offsets and the 
instrumentation is stable, the phase residuals for all scans (after subtracting the bin phase 
offsets) will be near zero, and changes in the phase variances caused by data drop-outs 
will have a minimal effect. The arbitrary line Fit chooses is one with slope equal to the 
average BSA delay for the band and has zero phase at the band reference frequency, cjo- 
As discussed below, this choice minimizes an error source in the final, charged particle 
calibrated delay. 

Obtaining Bin Reference Phases 

Before Fit can determine the bin phase offsets described above, it must first obtain 
the time-averaged phase points shown in Figure 5.1% This is done by fitting the phases for 
each bin in the band to the model phase d>"‘ pc , given by 

= ) + 4>T> l 5 - 1 ! 

where bin i has frequency a time point j has time tj , f mpc is the model phase delay rate 
for the band, and <f > is a reference phase for bin k (the points shown in Figure 5.1). The 
reference phases and the phase delay rate can be obtained by minimizing the chi-square 

*' = E i -2 I 5 - 2 ! 

ij 


5-3 



as a function of these parameters, where <f>ij is the phase of the data phasors and o\ 

J Yij 

is given by [2.15cj. Setting the partial derivatives of [5.2] to zero leads to the following 
equations, valid for all k : 


E ^ = fmp ' E + E [5.3] 

t t t 

D k = T mpc u k T k + 4i^S kt 

where the summations 


«»?;r 

J <Pij 

Y' (ty - to) 

‘ " 4- a* 

J Vij 

(*i~t o ) 2 


X ~^a\ 

J Vij 


5.41 


*-E 


J Yi} 


have been used. Solving equations [5.3] for <£g, pc and f mpc leads to 


f mpc 



4 


rape 

Oik 


Dk - 


[5.5] 


The expressions for £>, and in [5.4] involve using the phase of the data phasors. If the 
phase traverses cycles in time, as shown in Figure 7.1, D { and E { will have erroneous values 
because relative cycle ambiguities axe not assigned (the phase traversing several cycles 
during the scan is typical behaviour because the correlator model is usually not accurate 
enough to approximately stop the phase over a several minute scan). This problem is 
discussed and solved in detail in section 7. The solution there as here is to obtain accurate 
parameter estimates using Fourier transforms over time for each bin and counterrotating 
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these effects from the data. In this case, the fringe frequency, u/{, and phase, <f>\ x at time t 0 
are found from Fourier transforms for each bin. The data phasors are then counterrotated 

so that 

4>a = 4>a - <4 (*; - to) - <f>l M 


where 4> x j is the phase of the original data phasor and <j> x j is the counterrotated phase. 
If wj and 4>\ are determined with sufficient accuracy, the new average phase and phase 
delay rate will be about zero. If the system noise and unmodeled phase fluctuations are 
small, phase tracking (discussed in sections 12 and 13) is not necessary, and the following 
summations will be free of cycle ambiguity effects. 






fijj to) 
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Substituting these expressions and [5.6] into [5.4] gives the following relations 


D { = D { 4- wjTv + 4% [5.8] 

E x = E x + J X U X + 4>\T X . 

Finally, these are substituted into equations [5.5] to give 

£ u x {Ei — T x Dil S x 4 - u> x Zi ) 

•mpc __ _t 

t 

t> k + k - 

*ok = s; + 

where 

Z x = Ui - T?/Si. [5.11] 

This result depends on the two assumptions noted above: 1) the manual phase calibration 
(MPC) scan should have a high SNR relative to other scans in the experiment, and 2) 


[5.9] 

[5,10) 
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the phase residuals should not show deviations large enough to require phase tracking 
(this requirement could be eliminated if phase tracking were used to assign relative cycle 
ambiguities so the sums in [5.7] use the correct phase; Fit does not do this at present as 
this option has not been needed). The bin reference phases given by [5.10] correspond to 
those shown in Figure 5.1. The differences between these reference phases and a line with 
slope equal to the average BSA delay, 7" BSA , and with zero phase at wo are the bin phase 
offsets for bin k, <f>* k : 

$ k = K ~ r BSA (w* — Wo). [5.12] 

The BSA delay for each channel is calculated by performing a Fourier transform on all 
bin reference phases in that channel. This method gives the best estimate of BSA delay 
regardless of cycle ambiguities which are difficult to determine otherwise because of the 
large instrumental phase non-linearities across the 2 MHz bandpass. These channel BSA 
delays are then averaged over all channels in each band. 

When requested, Fit writes the bin phase offsets to an external file. These offset phases 
can be read in later and applied during the analysis of each scan in the experiment. This 
results in all phasors in a given bin being counterrotated by the same offset, for all scans. 
A possible improvement to the procedure given above would be to average the phase offsets 

from a number of scans; this would further minimize residual phases over the duration of 
the experiment. 

Note that since each observation band is processed separately, each has a separate 
arbitrary clock offset. This would pass through the charged particle calibration (section 
17) as a constant clock-like term if the delay reference band frequencies, given by [17.29], 
were constant for all scans; this is not the case in practice. The delay error from this 
effect is equal to the change in the delay weights, given by [17.6], over the scans in the 
experiment, times the difference in the S and X-band clock errors. The delay weights 
for a several hundred scan experiment are usually all within about 5 x 10“ 4 . As long as 
the difference in S and X-band clock errors is less than 20 nanoseconds, the delay error 
caused by this effect will be below about 10 picoseconds. In choosing the average BSA 
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delay as the line used to obtain the bin phase offsets, Fit attempts to minimize this error 
by minimizing the clock error. It should be noted, however, at this time, neither the 
instrumental ly induced clock differences nor the clock differences induced by using the 
average BSA delay as discussed in this section have been measured. Thus, the size of 
this effect has not yet been adequately determined (this error term is shown explicitly in 

equation [17.32a]). 

Manual Phase Calibration Consider Analysis: 

No Media Effects 

This subsection is devoted to obtaining expressions for the data phases after subtracting 
the bin phase offsets discussed above. The analysis presented here omits media effects but 
is expanded to include charged particle and tropospheric effects in the next subsection. 
The model phase, <£* m , is assumed to be given by 


tl” 1 = <t > o + r(w, — wo) + TWi(tj - to) + <f>l, [5.13] 

where <f>o is the residual (about the correlator model) phase, r is the residual delay, t is 
residual phase delay rate, f is an instrumental phase offset for each bin, assumed 
constant with time, j runs over all time points, tj , and i runs over all bin frequencies, 
Ui. The first three terms represent the phase without instrumental effects (see [6.2] and 
[14.1]). Substituting this phase into the expressions for D{ and Ei in [5.4] gives 


Di = [<f> o + r(wi — wo) + <^]S t ' + TU{Ti f 5 - 14 ] 

E{ = [4>o + r(wi - w 0 ) + 4>i\Ti + rWjI/j. 

Substituting these expressions into [5.5] gives 

r mpc = f [5.15] 

<t>T = <i>0 + r(w fc - wo) + 4>l [5.16] 
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Equation [5.15] indicates this method of obtaining the phase delay rate estimate is un- 
biased, again, assuming no media effects. Equation [5.16] is the expression for the Jfcth 
bin reference frequency. The bin phase offset for bin k, <f>* k , is obtained by subtracting 
^BSA(wj t — wo) from the reference phase (substitute [5.16] into [5.12]). Thus, the bin phase 
offset for bin k is given by 


= o + (r - TbsaKw* - w 0 ) + <f>{. [5.17] 

These are the phases which are subtracted from each time point in a given bin, for all 
scans in an experiment. If these MPC phases were applied to the same scan they were 
obtained from, the form of the resulting phases would be 

tij = <f>ij ~ <t>i [5.18] 

— r BSA (w t - — Wo) + TW t (iy — to), 

where the prime here denotes phases after MPC. This shows that applying MPC on the 
same scan from which the phase offsets were obtained leads to phases where the instrumen- 
tal phase offsets, <f>{, are removed, the band residual phase, <f> 0 , is set to zero, the delay is 
changed to the average BSA delay and the phase delay rate is unaffected (see the consider 
analysis at the end of section 14 for a proof of this statement). For other scans in the 
experiment, the calibrated phases are 

- <t>ij — $0 - i T * - r BSA)( w t - W 0 ) - (j>[ [5.19] 

= [<t> 0 — <j> o) + (r — T* + Tb SA )(w,- — Wo) + TU>i(tj — to), 

where the * superscript denotes quantities from the MPC scan. This shows the resulting 
band residual phase will be offset by <f>*, the delay will be offset by the constant clock-like 
term equal to r* - r^ SA , and the resulting phase delay rate will be unaffected. In sum- 
mary, if media effects are neglected, manual phase calibration eliminates the instrumental 

phase offsets and induces a constant clock-like term and phase offset for all scans in the 
experiment. 
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Manual Phase Calibration Consider Analysis: 

Dispersive and Tropospheric Effects 

The following subsection expands the consider analysis to include dispersive and tro- 
pospheric effects. This material can be skipped on first reading; much of it is related to 
material found in section 17 on charged particle calibration. To avoid repetition, a discus- 
sion of dispersive effects on phase, the meaning of fcy below, the sign convention and other 
supporting material can be found in section 17. 

To include dispersive and tropospheric media effects into the consider analysis, equation 
[5.13] must expanded: 

</>£ = 4>o + T{ui - W 0 ) 4- - to) + <j>\ - fcyM + T j U i’ l 5 - 2 °] 

where Jfey is proportional to the differential total electron content in the Block II sign 
convention and rj is the differential tropospheric delay at time ty. Expanding kj and 7y 
about a linear model gives 


Jfey = k + (ty - t 0 )k + fcy l 5 - 21 ] 

r j = + (ty - to)f l + ry r , 

where all quantities are in the Block II sign convention. Substituting these back into [5.20) 
gives the model phase including dispersive, tropospheric and instrumental effects: 


4>j - = (<£o + r'wo) + {t + f‘)(wi - w 0 ) + {t + T)oJi{tj - t 0 ) 


[5.22] 


, k k(tj to) 

+ 4 ' i ~7 i ^ 


h w,ry r . 

w.- 1 


Note the linear component of the tropospheric delay enters the first three groups of terms 
as expected; the phase and group delays have an additional and the phase delay rate 
has an additional t* (see the consider analysis at the end of section 14). This equation can 
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be substituted into the expressions for Z? t and E { in [5.4] to give analogous expressions to 
[5.14]: 


A = [(*> + + (r 4- - w 0 ) + + (r + [5.23] 

kS t / kTi/tjJi — Fj/u>i + 

A = [(<£o + r^o) + (r + f‘)(w t - - wo) + ^]T t - + (f + r f )w. A 
kTijwi — kUi/u>i — G</w f + ujjNj, 


where 


k r 


r* r 


<71 




<71 

J 0»J 


[5.24 


«<■£ 


*o) 




rtj 


rxj 


Substituting equations [5.23] into [5.5] gives expressions analogous to [5.15] and [5.16]: 


r m pc = f + f* - k ~ FiTi/Sj) + — MiTjj Si) 


32iufZi 


Hi wfz t 


5>?*i 


*5T — 4>o + tuq + (r + - w 0 ) + <j>{ Ic 


Tk 

“k v k S k 


1 - 


(jj 


k Hi Z\ 


Hi ™}Zi\ 


_ jjfc _j_ HiiHi ~ 

“JfcSjfc 5* ^wfz, 

+ ,.,*** w*T* “, 2 (A - M(Tt/St) 

k s k Sk Hi"fZi 

The bin phase offsets can be found as above to give 


[5.25] 


r k = 


WS + f, '"5) + ( r * + - «•;«)(»>* -«{) + #{- — 


+ ( £ 3A + <5* - <6*), 


[5.26] 
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where wj is the MFC scan’s reference frequency, the following time definitions, 



3U, i l 



[5.27] 


2 ‘ " 

,, n e,^ - MZfi) 

3 ‘”s t 


.. = Mk 

hk ~ S t 

. n r^UN.-MMSi) 

6 t ~s k EjUijZi 


are used, and the superscript * has been inserted to explicitly show quantities from the 
MPC scan ( f tt , the analogous expression to t ik for the tropospheric delay, is zero and will 
not be used elsewhere in this report). 

Note that t^ sa is corrupted by dispersive effects in the MPC scan. Dispersive effects 
corrupt the BSA delays for each channel, which are averaged to obtain r^. Introducing 
a small, constant correction delay, r c * or , the BSA delay can be written 


'BSA 


= r. 


+ fc*K 2 + r* 


BSA 


cor > 


[5.28] 


where u* (including the band subscript) is given by [17.29] for the MPC scan and r* SA 
is the band averaged BSA delay if no dispersive effects were present. The k*/u? term 
accounts for almost all of the dispersive effect; r‘ or is only introduced to make [5.28] a strict 
equality. The bin phase offsets can now be written 


■K = « + + K. + f< * - )("* - "») + +1 - 1 


k* 


U k 


[5.29] 


7*tU 


— k — — — + + *lk ~~ *6 *)’ 

u k U k 


where 


* * 

'cl 


— T r BSA r cor 


[5.30] 


both represents the constant clock-like error discussed above and is used to absorb the 
small correction term in [5.28]. 
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If MPC is applied to any scan in general, the phase after calibration can be written 
(see [5.19] for comparison): 


<t>ij = <t>o ~ 4>l + wo (r 4 - r 4 *) - Au; 0 (r* { - ] 


5.311 


+ (r - r c * + f 4 - f 4 * 


k* 


+ ~ w °) + ~ *o) 

tl . . 


- k ~ir - lih ~ + v % + -e*2f - + - <«> - - * + . 


where 


Aw 0 = w 0 - u/q. [5.32] 

The first two lines of terms show the effects on the phase, group delay and phase delay 
rate, the next six terms describe dispersive effects and the last term describes the effect 
from the non-linear component of the tropospheric delay. Note the mixture of effects from 
both the scan in question and the MPC scan. A consider analysis in section 17 will show 
the effect each term has on the phase fit observables, the group delay, phase delay rate 
and phase delay, for each band. 
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Section 6 


A Priori Parameter Estimation 

This section, along with sections 7 through 11, is concerned with source detection, 
defined in section 10, and obtaining an a priori estimate of the residual BWS delay, phase 
delay rate and phase for each observation band (the band phase observable is usually not 
important except when phase connection is used to derive the phase delay, as given in 
[18.1]). These a priori estimates are obtained to effectively improve the correlator model; 
the model phase derived from the a priori observables is removed from the data to simplify 
further processing. The model phasor used for source detection, P- , is given by 

= A b e< (6.1] 

where i runs over all bin frequencies w, in the band, j runs over time points tj, A b is the 
band amplitude, <£“■ is the model phase given by 

<t>“j = 4>o + r a K - wo) + - to), l 6 - 2 ) 

T a is the a priori residual BWS delay, T a is the a priori residual phase delay rate and 4>q 
is the a priori residual phase. The term “residual” here is relative to the correlator model 
which has already been subtracted from the data phasors by the Block II correlator. The 
parameters A b , 4 >q, r a and f a are to be estimated by coherently integrating the data as 
described in the next several sections. After determining these preliminary values, the 
phase of each data phasor is counterrotated by the model given by [6.2]. If the signal 
strength is great enough, the resulting counterrotated phases can be tracked as discussed 
in sections 12 and 13 and fit using the linear least-squares analysis described in section 
14. If the signal strength is too low, the counterrotated phasors are fit using the linear 
least-squares analysis discussed in section 15. 

The band amplitude, given by A b in [6.1], is obtained from the a priori estimation 
procedure but is not the final amplitude Fit passes to Modest. Equation [14.17] gives the 
estimated band amplitude sent to Modest (see reference 9 for more details). 
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It is important to note that although the model given by equation [6.1] is useful for 
source detection and obtaining the a priori estimate discussed in the next few sections, it 
is insensitive to phase changes of an integral number of cycles. A loss of coherence will 
result if the residual phase about this model fluctuates over several cycles. In practice, 
the coherent integration implemented with [6.1] will detect moderately weak signals even 
if cycle fluctuations on the order of several cycles are present but will result in inaccurate 
a priori parameter estimates, especially for the phase delay rate. This is a common prob- 
lem at S-band (and even worse at lower frequencies) where ionospheric charged particle 
fluctuations commonly induce phase changes on the order of a cycle per minute. Larger 
fluctuations are seen in observations less than about 15 degrees from the sun due to solar 
plasma fluctuations. If the detection SNR is high enough, phase tracking with a sufficiently 
high sampling rate can improve the parameter estimation. This is equivalent to an im- 
proved model because phase tracking, defined in section 13, is sensitive to phase changes 
of an integral number of cycles. If the SNR is high enough for detection but too low to 
allow phase tracking, no verification that the phase has not traversed cycles can be made, 
making the delay and phase delay rate estimation more inaccurate than the formal error 
would indicate. 

In summary, the model given by equations [6.1] and [6.2] is used for the a priori 
estimate because of its ability to detect weak sources and because it allows the fast Fourier 
transforms discussed in the following sections. If the SNR is high enough, phase tracking 
can be used to improve the model and thus the parameter estimates. If the SNR is not 
high enough to allow phase tracking and there is some possibility of phase fluctuations 
from any source being on the order of one cycle or more, the formal errors on the final 
observables will be underestimated. 
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Section 7 


Time Fourier Transform to Obtain 
Residual Fringe Frequency 

The Block II correlator passes to the user the complex correlation coefficient, given by 
[2.13], at each dump interval and each data bin, typically every two seconds and four per 
2 MHz channel, respectively. These data are fit to the model given by [6.1] and [6.2] to 
extract the a priori BWS delay (r tt ), the phase delay rate (f a ) and phase (<*§). This would 
be a simple least-squares fit in the absence of cycle discontinuities in the phase. Figure 7.1 
shows the phase of the complex correlation coefficient as a function of time for a single 
data bin. The phase can be seen to traverse several cycles over the scan, where the slope of 
the line passing through these points is the residual fringe rate. This behavior is expected 
because the correlator model is generally not accurate enough to stop such a small fringe 
rate: about 0.025 Hz in this case. A naive linear least-squares fit to the phases would not 
produce the solid line shown in figure 7.1, but rather a line with approximately zero slope 
along zero phase. Thus, cycle ambiguities must be resolved before least-squares fitting. 

There are simple algorithms which allow cycle ambiguity resolution. For example, 
moving each point in the time series to within a half cycle of the previous point would line 
up the phases as shown in Figure 7.2. This algorithm has the disadvantage that the signal 
must have a sufficiently high SNR and low residual fringe rate that there is little chance of 
a time point being truly more than half a cycle away from the previous point; otherwise 
the algorithm will erroneously move the next point and thus all later points in the time 
series. When the SNR is moderate to low, or the residual fringe rate is large, a search 
algorithm is the only robust approach to finding the residual fringe rate and phase. One 
could search for the best solution by fitting the phase time series to a large range of slopes 
(residual fringe rates) and intercepts (phases), while handling the cycle changes explicitly. 
This is essentially the approach taken here except, without the refinements given below, 
it is computationally too expensive. 

Imagine subtracting a linear phase model of the form (f>^ m {t) — ft t from the residual 
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Ftg. 7.1 Residual phase as a function of time for one bin. Because the 
phase traverses several cycles, a simple linear least-squares fit to 
these data cannot be performed. The solid line shows the best 
fit if cycle changes are handled correctly. 



Fig. 7.2 Residual phase as a function of time after cycle ambiguities have 
been resolved. The solid line shows the linear least-squares fit 
to the data. 


phase shown in figure 7.1, where n is a trial fringe rate. If 0 were close to the true residual 
fringe rate, the phases would line up almost horizontally, as shown in figure 7.3. A search 
for the n which results in the best horizontal alignment is equivalent to searching for the 
best residual fringe rate. The concept of subtracting a trial phase model and looking for 
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Time (Seconds) 

Fig. 7.3 Residual phase after subtracting a trial linear model. The solid 
line shows the linear least-squares fit to the residual phases. 

a constant phase can be extended to the complex correlation phasors. Here, the phasors 
are counterrotated by multiplying by e~ int . If fl is close to the residual fringe rate, the 
phasors will, on average, line up with the same phase. One can test for this condition by 
calculating the amplitude of the phasor sum, which can be expressed as 

A* = , M 

where p fc (t y ) is the phasor at the jth time point, given by [2.13]. Only when the phasors are 
lined up will they sum coherently, resulting in a large relative amplitude. This phasor test 
has the advantage that it does not fail when the phase is very close to ±180 degrees, where 
the average phase can be calculated to be near zero (if about half the phases are around 
+180 and half around -180 degrees). The amplitude can be calculated for a number of 
trial frequencies, H*, as 

4= i>(*, i 7 - 21 

J=1 

so that the maximum amplitude as a function of k corresponds to the best trial residual 
fringe frequency. 
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Equation [7.2] is quite similar to that of the discrete Fourier Transform 




n-l 


3=0 







3 = 1 


[7.3] 


where k = 0, 1, . . . , n - 1 (see [B.l] and [B.2] for FFT definitions). Equating the exponents 
in [7.2] and [7.3] gives the effective frequency tested as a function of k. Performing this sum 
explicitly is still computationally expensive; the number of computations grows as n 2 . If n 
is a power of two or can be factored into powers of small primes, a fast Fourier Transform 
algorithm (FFT) can be used to greatly cut the computation time. If n is a power of two, 
the number of computations grows approximately as n log n. The FFT can be improved 
for this application by defining additional phasors with zero amplitude beyond the actual 
number of data phasors. This padding of the FFT algorithm with zero-amplitude phasors 
results in a finer resolution in the transform space (residual fringe frequency), but does 
not change in any way the result for a given value of residual fringe frequency; it only 
decreases the spacing between trial values (see appendix B for details). 

Fit performs an FFT on each data bin using n chosen to be a power of two between two 
and four times the number of data phasors. Such an n is greater than necessary to uniquely 
determine the maximum-amplitude phasor and obtain the residual fringe frequency, but 
this higher density of points improves the rescaling and interpolation discussed in the next 
section. 
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Section 8 


Rescaling Residual Fringe Frequency 
to Residual Phase Rate 

Transforming each bin with an FFT over time results in a new set of phasors. The 
new phasor with the greatest amplitude corresponds to the best estimate of the residual 
fringe frequency for that bin, given the discrete sampling of the FFT; the true maximum 
may lie between two transformed phasors. Coherently summing the data across all bins 
in a frequency band both increases the SNR (necessary for weak sources) and allows the 
determination of the a priori residual BWS delay, r a . Such a sum involves a second FFT 
across frequency, discussed in the next section, but should not be performed at this stage 
of the analysis due to the residual fringe frequency dependence on bin frequency. The 
physical observable desired is the residual phase delay rate, f, and is equal to the residual 
fringe frequency, w/, divided by bin frequency; r = u/M, where w t - is the observation 
frequency of the ith bin. Because r is constant over bins, the maximum phasor for each 
bin will occur at a slightly different fringe frequency. Thus, a coherent sum over bins 
should not be done until the peak amplitudes are aligned to avoid amplitude smearing. 
This effect is shown in the upper plot in figure 8.1. 

If wo is the band reference frequency, a scale factor V’i = o can be defined so that 

(iff = TU>i — fljjQlpi = [8.1] 

where <[> = ojot is the phase rate. The bin frequency dependence of the transformed phasors 
can be removed by contracting the scale of residual fringe frequency by ipi for each bin, 
resulting in a new phase rate scale. Computationally, nothing is done except to reinterpret 
the residual fringe frequency associated with each transformed phasor. The middle plot 
in figure 8.1 shows the effect of rescaling; it has changed the spacing between points in 
different bins. The final step before the coherent frequency sum is to calculate phasors 
corresponding to a common set of phase rates for all bins in the band. This is done by 
interpolating the phasors already in hand to these common phase rates. 
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Fig. 8.1 Schematic representation of fringe frequency rescaling. The up- 
per plot shows an exaggerated example of the phasor amplitude 
after the time FFT for four bins. The fringe frequency depen- 
dence on bin frequency causes a shift in the maximum amplitude. 
The middle plot shows the effect of rescaling the horizontal axis 
to phase rate. The lower plot shows the interpolated values 
along a regular grid so that an FFT over frequency is possi- 
ble. The density of points reflects that used by Fit; the top 
and middle plots have a higher density to allow a better pha- 
sor interpolation, but the density is reduced in the lower plot 
to satisfy computer memory constraints. The plots on the right 
show more clearly the relative spacing of points in the plots on 
the left. 
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The phasor interpolation results in a phasor with amplitude equal to the weighted 
geometric average and with a phase equal to the weighted arithmetic average of the two 
interpolating phasors. This can be calculated quickly using the following 

| 8 - 2 | 

where p is the interpolated phasor, p\ and p 2 are the interpolating phasors and h is the 
difference in phase rate between p and p\ divided by the phase rate difference between P 2 
and p\. The actual calculation counterrotates both phasors by the phase of the nearest 
interpolating phasor before the multiplication, and rotates the result back the same amount 
to avoid the branch cut along the negative real axis when doing the complex exponentiation. 

As noted in the previous section, to improve the phasor interpolation accuracy, the 
density of fringe frequency points across each bin is about twice that needed to find the 
maximum. The program Fit chooses the common set of phase rate points to have a 
density half that found in a hypothetical bin having the band reference frequency, ujq. 
This reduction in the density of points is necessary to satisfy computer memory constraints, 
but is dense enough to locate the maximum unambiguously. The final plot in figure 8.1 
schematically shows the interpolation and density reduction. The interpolated phasors are 
now ready to be coherently summed over frequency. 
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Section 9 


Frequency FFT to Obtain 
Residual BWS Delay 


The cross correlation phasors, after having each bin’s time series transformed to resid- 
ual fringe frequency and rescaled to residual phase rate, are ready to be combined to 
increase the SNR and derive the a priori residual BWS delay. This combination is done 
with an FFT over frequency. 

A possible procedure at this point, though not the one used by Fit (the first, private 
version of Fit did do this), is given to motivate the overall procedure. After the time 
FFT and rescaling are performed on a high SNR observation, the maximum-amplitude 
phasor in each bin can be used to give the best estimate of the residual phase rate, and 
thus phase delay rate, for that bin; a fit to interpolate between the discrete FFT phasors 
is used to obtain the best estimate. A weighted average of these rate estimates from all 
bins can give an overall estimate of the residual phase delay rate for the band. This can 
then be multiplied by each bin’s observation frequency to give the best estimate of the 
residual fringe frequency for each bin. The phase of the original time series of phasors 
(before the time FFT) for each bin can be fit to a line with slope equal to this calculated 
residual fringe frequency, with the fit determining the phase intercept at the scan reference 
time, to (this assumes the SNR is great enough that phase tracking in time can be done 
easily). Figure 9.1 shows an example of these phase intercepts plotted as a function of bin 
frequency (this procedure is very similar to that discussed in section 5 on manual phase 
calibration). The slope of the line connecting the points is the residual BWS delay, r°, 
and the phase intercept at the reference frequency is the residual phase, d>Q, in equation 
[6.2]. The slope of the line connecting the bins for a single channel is the BSA delay for 
that channel. In analogy with the time sequence of points shown in figure 7.1, the cycle 
changes in figure 9.1 make fitting a line difficult. The only robust method of determining 
the slope is by fitting with trial slopes and searching for the best fit. Just as in the time 
series case, an FFT is needed to perform the search in an efficient manner. 


9-1 



o 



F tg. 9.1 The fit phase at the reference time as a function of bin frequency. 

Because the phase traverses several cycles, a simple linear least- 
squares fit to these data cannot be performed. The solid line 
shows the best fit if cycle changes are handled correctly. 

The analysis presented in the last paragraph has the disadvantage that low SNR ob- 
servations will preclude the time series fit to the original data. Performing the frequency 
FFT directly on the time-transformed phasors avoids this problem. The resulting two- 
dimensional (2-D) FFT is equivalent to a 2-D search for best residual BWS delay and 
phase delay rate estimates; the amplitude of a phasor in the 2-D transformed space indi- 
cates how well that trial delay and rate fit the data. 

The FFT over frequency has one complication not found in the time series FFT: the 
points to be transformed using the fast algorithm must occur at integer intervals but 
the bin frequencies are often selected by the experimenter to avoid even frequency gaps. 
This problem can be approximately solved by noting that the lowest frequency channel in 
the band has n bins which can be assigned to the first n integers. This sets the slope and 
offset needed to calculate the integers corresponding to the other bins in the band. If other 
bins do not lie exactly on an integer, the nearest one is chosen. This is usually a good 
approximation because the channel LO frequencies are typically separated by multiples of 
0.5 MHz with small offsets on the order of 0.01 MHz. This approximation is later dropped 
by calculating the Fourier Transform explicitly as discussed in the next section. Integers 
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not associated with actual data points use zero-amplitude phasors, and additional zero- 
amplitude phasors are added to the end of the sequence to at least double its length and 
to make the total number of points a power of two. 

After the phasors are Fourier transformed over frequency, the original two-dimensional 
time/bin frequency data will be transformed to a residual phase delay rate/residual BWS 
delay plane. A search for the maximum-amplitude phasor will correspond approximately to 
the best estimate of residual phase delay rate and residual BWS delay; the best estimate 
corresponds to a point near the maximum phasor due to the discrete sampling of the 
FFTs. The phase of the maximum-amplitude phasor corresponds to the residual phase, 
4 > q . The two-dimensional grid of points is usually dense enough to unambiguously find the 
nearest phasor to the true peak amplitude among the discrete transformed phasors, but is 
not dense enough to obtain sufficiently accurate estimates of the residual observables. A 
higher density of points could be obtained by padding the FFTs with more zero-amplitude 
phasors, but this option is eventually restricted by memory limitations. The solution used 
by Fit is to perform the two-dimensional Fourier sums explicitly to obtain the transformed 
phasor at any point; an iterative search is used to locate the true maximum-amplitude 
phasor (see Appendix B). The next section covers this technique in detail. 
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Section 10 


Exact Fourier Transform to Obtain 
Residual Phase Delay Rate and BWS Delay 


After applying the time FFT, rescaling, and frequency FFT, the bin/time phasors are 
transformed to residual phase delay rate/BWS delay phasors. The phasor in this latter 
space having the maximum amplitude corresponds approximately to the best estimate of 
the residual phase delay rate and BWS delay. The phase of this phasor corresponds ap- 
proximately to the best estimate of the residual phase. Because memory and computation 
constraints limit the size of the residual phase delay rate/BWS delay grid, an interpolation 
is needed to find the true maximum-amplitude phasor. One approach is to use the phasors 
already calculated to fit the data with the expected functional form near the peak. This 
function is approximately sinx/x in the rate direction and the delay resolution function 
in the delay direction. This approach is rejected because 1) the complicated functional 
form makes the fit difficult and 2) when large temporal phase fluctuations are present, 
this is not an accurate model anyway. The method used by Fit is to calculate the trans- 
formed phasors explicitly using a slow Fourier transform. This allows the calculation of 
the transformed phasor at any point in the residual phase delay rate/BWS delay plane. 
An iterative search for the maximum amplitude is then performed. This approach has the 
added advantage that the approximations made to obtain the phasors in hand, namely the 
phasor interpolation used in the rescaling operation and the rounding of the bin frequen- 
cies needed for the frequency FFT, can be dropped in the explicit calculation. In addition, 
the iterative search utilizes a fairly straightforward and robust algorithm. The primary 
disadvantage of this method is the calculational cost of performing an explicit FT. 

The explicit FT is designed to reproduce exactly those phasors already calculated with 
the FFTs, assuming the two approximations mentioned above have a negligible effect. The 
detailed steps needed to derive the form of the explicit FT can be found in appendix B, 
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but the final equation can be written as (see [B.9]) 


f„ = V ex P 2rif (t ~ 1)(2 ‘~ n,/2 ~ 1) ^- + 

nTfr? L n t W 0 uj b n f J 


[ 10.11 


where pu is the data phasor for the Arth time point and Ith. bin, T tJ is the transformed 
phasor for the tth phase delay rate and jth BWS delay point, n b and n r are the number 
of bins and data time points respectively, nj and nt are the number of BWS delay and 
phase delay rate points in the respective FFT and are each a power of two, and wo, uq 
and wjfe are the band reference frequency, the first bin’s frequency and fcth bin’s frequency, 
respectively. Comparing the exponent phase in equation [10.1] with the model phase in 
equation [6.2] leads to the following identifications (see the last subsection in appendix B 
where these equations are derived in more detail): 
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where t c is the interval between time points, usually two seconds, <j> m&x is the phase of the 
maximum-amplitude phasor, to is the scan reference time, t\ is the time tag of the first 
point in the time series, u) b is the bandwidth of a single bin and 
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[10.3] 


The variables x, and Xj correspond to t and j in equation [10.1] but can take real values 
as opposed to integers (necessary for the interpolation). 
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The first step in the search for the maximum-amplitude phasor is to recalculate the 
transformed phasors with the exact FT for a three by three grid centered on the maximum 
phasor found with the FFTs. This is done to ensure that the maximum does not change 
location when dropping the two approximations noted above. In some rare cases the peak 
will move, usually by just one point (if it does, a new three by three grid is calculated). 
After locating the maximum amplitude in the three by three grid, the center phasor and 
the nearest four adjacent phasors are used to perform a two-dimensional parabolic fit. The 
amplitude of the five phasors is fit to the function 


A m = oii^ 4- 0 , 2 Xj + a3Xj + a\Xj 4- as 


( 10.41 


where the coefficients ai through a 5 can be determined from the center and nearest four 
phasor amplitudes. This paraboloid has a maximum at the point 
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where the superscript c refers to the location of the central phasor, Ax* and A Xj are the 
spacings between the phasors in the X{ and xy directions, and the terms are the five 
phasor amplitudes arranged in the two-dimensional X{ and xy plane as follows 

A m 

^ 0 + 

A m Am Am 
^ i _0 ^00 ^+0 ’ 

Am 
A 0 - 

The center of the paraboloid is taken to be the new maximum amplitude location and five 
new phasors are calculated: one centered on this point and four surrounding it as shown 
above. Typically, the center phasor has the greatest amplitude of these five. If this turns 
out not to be true, the phasor having the greatest amplitude is chosen to be the center 
and the four phasors adjacent to it are calculated. The spacing between the phasors in the 
grid is also reduced to converge on the maximum more quickly. This process is repeated a 
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number of times until the last iteration gives the final location of the maximum-amplitude 
phasor. 

If the amplitude of the final, maximum-amplitude phasor is significantly greater than 
the amplitude of phasors away from the peak region, the peak is assumed to be a result of 
the source signal; this defines source detection. If no source is present, or the source is too 
weak to detect, a maximum amplitude phasor will still be found but it will usually have 
an amplitude only about 3 times greater than the average amplitude of phasors away from 
the peak region; this ratio is usually at least 10 for well detected sources. 

The maximum-amplitude search just described is very CPU intensive. To reduce com- 
putation time, the number of iterations and the grid spacing reduction rate have been 
empirically optimized. At present, Fit reduces the grid spacing by a factor of four on each 
of four iterations. A larger reduction in grid spacing results in a greater probability that 
the center of the five new phasors does not have the greatest amplitude. This essentially 
adds at least one iteration to the nominal four. Even with this optimization, about 25% of 
Fit’s total processing time is used here. Shifting some of this computational load back to 
the FFTs by padding the time and frequency arrays with more zeros is less efficient. For 
example, doubling the number of points in the time and frequency FFTs will increase the 
computation there by about a factor of four, which at present, accounts for about 30% of 
Fit’s total computation time. 

Once the maximum-amplitude phasor is found as a function of x t and xy, and its phase 
<£m»x calculated, equations [10.2] are used to obtain the best a priori estimates of residual 
BWS delay, r«, residual phase delay rate, f*», and residual phase, No approximations 
have been made in deriving these quantities, and they are exact if the iterative search 
for the maximum phasor is done with sufficient accuracy and the model shown in [6.1] is 
appropriate. Typically, equation [6.1] is not the best model as discussed in section 6 and 
further analysis is necessary as described in sections 12-15. 
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Section 11 


Phase Counterrotation by the 
A Priori Model 

After deriving a priori estimates for the residual BWS delay, phase delay rate and 
phase, the phase of the data phasors is counterrotated by the model given in equation 
[6.2]. This is done to allow more robust phase tracking algorithms, discussed in the next 
section, and to permit an initial view of the phase residuals with Fit’s plotting software. 
If pij is the original data phasor for frequency u>i, at time tj , the counterrotated phasor, 
p c , can be written 

Pfi = |n i] 

and <t>*j is identical to [6.2]: 

<f> jj = 4>0 + T a {ui — Wo) + f a U{ (tj — to). [H-2] 


Parameter estimates based on these counterrotated data are now residuals of the a priori 
estimates, which are themselves residuals of the correlator model. The term residual should 
now refer implicitly to the estimated parameters which are added to the a priori estimates 
and the correlator model to obtain total observables. 
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Section 12 


Phasor Time Integration 

This section and the following three are concerned with obtaining the final estimates 
of residual BWS delay, phase delay rate and phase using least-squares analysis. Accurate 
estimates of these parameters have been determined and removed (counterrotated using 
[11.1] and [11.2]) from the data using the phasor model given by [6.1] and [6.2]. If the 
SNR is great enough, as defined below, the phasor model can be improved by tracking the 
phasor’s phase as a function of time (discussed in the next section) and fitting it to the 
model given by [6.2] alone (section 14). If the detection SNR is too low, a linear least- 
squares fit using the phasor model is performed to obtain the residual parameter means 
and covariance matrix, as discussed in section 15. The parameter means for this latter fit 
axe small, and would be zero in the limit that the iterative search described in section 10 
were done perfectly, because the model is identical in both cases. This section is concerned 
with determining which analysis path, phase tracking with the phase fit or the phasor fit 
alone, should be used. 

Phase tracking involves assigning an integer number of cycles to the phase of each 
data phasor in the observation band so that the relative phase difference in each bin’s time 
series and each time’s frequency series is handled correctly; an overall cycle ambiguity 
for the entire band is left unresolved (phase connection, discussed in section 18, involves 
determining the overall band cycle ambiguity). Some assumptions about the phase time 
variation must be made to correctly assign the integer cycles to each data phasor. The first 
assumption is that over a correlator integration interval the phase residuals fluctuate less 
than about one radian. Unmodeled fluctuations in the phase delay rate greater than about 
a radian per integration interval cause coherence loss, leave the relative phase ambiguity 
between time points difficult to assign and create a phase delay rate error greater than 
that predicted by the formal error. Tracking the phase with more frequent correlator inte- 
grations can help solve these problems but at the cost of lowering the SNR per integration. 
If the source SNR is too low, sampling the phase often enough to track the fluctuations 
will fail because the phase error on each phasor will be so large that the cycle ambiguity is 


12-1 



again arbitrary. The strategy used by Fit is for the user to obtain correlator integrations 
as often as possible: usually every two seconds with the present Block II configuration 
(with a much greater effort, the standard Block II sampling bandwidth of 0.5 Hz can be 
raised to 40 Hz). If the phase errors for these highest bandwidth phasors are too large, 
Fit time integrates the phasors in software until a sufficient fraction of them meet an SNR 
threshold; the small fraction of data points not making the SNR cutoff is deleted. This 
results in a phasor time series with the greatest sampling bandwidth, with each phasor 
having at least the threshold SNR. If the phase fluctuates by more than about a radian 
over one of these software integration intervals, the phase tracking may still give erroneous 
results; either the data are not sensitive enough to reliably track the phase or the SNR 
threshold is too high. 

The phasor SNR is defined to be R/o, where R is the center of the phasor’s amplitude 
probability distribution. This is estimated using A/ a (since R is an unknown) where A 
is the measured phasor amplitude and a is given by [2.9] or [2.19]. Figure 12.1 shows 
the phase probability density, given by [2.15a], as a function of phase for different SNR 
levels. As the SNR level drops, the tails of the distribution from different cycles begin to 
merge. This results in a much higher probability for an outlying phasor to have its cycle 
ambiguity resolved incorrectly, even when the surrounding phasors in the time series have 
small unmodeled phase fluctuations. The plot indicates the SNR threshold should be about 
two; a higher level lowers Fit’s software sampling bandwidth, while a lower level increases 
the probability of resolving the relative cycle ambiguities incorrectly. Empirically, the most 

efficient value of the SNR threshold changes from scan to scan but is often between 1.5 
and 2.0. 

Phasor Integration 

The phasors are integrated in software using a variation of equations [2.13] and [2.14]: 

- _ ( Cl S[ \ 

1 '7V//4-I//2’ JVj/4- J//2/’ f 12 ' 1 ! 

where p/ is the resulting integrated phasor, S Iy C,, and // are the like-sign sine, like-sign 
cosine, and invalid counts and Nj is the total number of bits, summed over all phasors in 
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Probability Density 



Cycles 

Fig. 12.1 Phase probability distribution for various SNR values; see equa- 
tion [2.15a] and [A. 14]. The probability distribution for high 
SNR levels is nearly Gaussian and remains clearly separated 
between cycles. For SNR levels below about 2.0, the distribu- 
tion flattens and the contributions from adjacent cycles become 
merged. 

the integration interval. Similarly, the error on each phasor component is given by 



Vf ~ 2 Bj_ 
N l V,'‘ 


[ 12 . 2 ] 


where Vj is the integrated valid counts, Vj = iVj — 2//, and i?j is the integrated blanking 
counts. 
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Fit incrementally increases the integration interval and calculates the fraction of data 
phasors having the threshold SNR. The smallest integration interval in which a specified 
fraction (usually about 0.90-0.95) of integrated phasors exceeds the SNR threshold is used. 
If the integration interval exceeds | of the scan’s duration or a user-specified time interval 
(usually 1 minute), whichever is smaller, without the specified fraction of phasors having 
the SNR threshold, Fit does no phase tracking and fits with the phasor fit. Otherwise, 
the sequence of integrated phasors is phase tracked as discussed in the next section and fit 
with the phase model. For long scans, the upper limit on the time integration should be 
less than \ the scan duration because the chance of integrating over cycle slips increases 
with the integration interval. 

If tone data are used to calibrate the cross-correlation phasors, as described in section 
4, the tone accumulator counts must also be summed. The summed tone counts are used 
to calculate the tone phasors using an expression analogous to [4.1]. The summed cross- 
correlation phasors are then calibrated using equations [4.2] and [4.3]. 

Shift in A Priori Estimates 

Before further analysis, the final values of t 0 and U7 0 are computed. The phasor sums 
and the elimination of summed phasors not having the SNR threshold have changed both 
the phasor time tags, which are always centered in the integration interval, and the phase 
variances, a\... Equations [3.1] and [3.2] are re-evaluated at this point in the analysis to 
give the final t 0 and w 0 values. This re-evaluation induces an implicit change in the a priori 
estimates; to keep the counterrotation phase <£“■, given by [11.2], unchanged, shifts in the 
a priori estimates are necessary to compensate for changes in t 0 and w 0 . If primes denote 
new values, the counterrotation equations become 


+ r °( w t ~ Wo) + f a u>i(tj — to) 
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Setting <£“ equal to <£?■' and equating terms with the same t and j dependence leads to 


<f>o = <£o + ( w o — ^o)r a + (io — t o)(x>' 0 T a [12.4] 

T a> = r a + (t' Q - t 0 )r a 


Fit uses these expressions to adjust the a priori estimates; further references to a priori 
estimates imply these modified (primed) values implicitly. 
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Section 13 


Phase Tracking to Improve the Model 


After integrating the phasors in software so that each resulting phasor has at least the 
threshold SNR at the highest sampling bandwidth, the cycle ambiguities for each phasor 
are assigned. Figure 13.1 shows the residual phase as a function of time for one of several S- 
band channels after removing the a priori phase estimate. The phase clearly traverses two 
full cycles near the center of the scan. Because the a priori model is insensitive to integer 
cycles, this solution was presumably obtained for its linear regions near the beginning and 
end of the scan. If the phasor model, [6.1], were used to fit the residuals, this would be very 
close to the best fit solution because the models are identical. This solution results in an 
incorrect phase delay rate; the true slope after the phase ambiguities have been resolved 
can be seen by eye to be about 0.014 Hz (2 cycles over 140 seconds), corresponding to 
a bias in the phase delay rate of about 6.1 x 10 -12 . The error in the delay direction is 
more difficult to calculate; the phases of the other channels in the band are necessary to 
completely analyze the BWS slope (Monte Carlo simulations using a random walk for the 
phase fluctuations indicate that fluctuations on the order of a cycle per minute result in 
a delay error on the order of 10 picoseconds). If the SNR had been too low to observe 
the cycle changes, or if the phase changed too fast for the correlator sampling to track 
(two seconds per point in this case) , the large bias in the rate and perhaps a small delay 
bias would go unseen, leading to an underestimate of the formal errors. A phase fit at 
this point (with no phase tracking) will definitely result in a biased delay and rate; the 
resulting rate would have the same problem as the phasor fit and the delay near a cycle 
crossing, for example at t = 95 seconds, is large and biased because the phase in the first 
three bins is near -b 180 degrees while the last bin’s phase is near —180 degrees. 

Resolving the phase ambiguities in figure 13.1 can be easily accomplished by eye, but 
is much more difficult to implement in software. Several empirical algorithms have been 
tried and rejected. The algorithm presently used by Fit performs reasonably well in most 
cases. The algorithm makes the assumption that the unmodeled phase fluctuations are 


13-1 



RESIDUAL BIN PHASE. CHANNEL 1 



Fig. 1S.1 Residual phase as a function of time for the four bins of a single 
channel. The scale for each bin is ±180 degrees. The phase 
traverses two cycles near the scan center. 

approximately equal for each channel in a band. This is not strictly true for dispersive 
effects where the phase change is proportional to the inverse observation frequency or for 
tropospheric delays where the phase fluctuation is directly proportional to the observation 
frequency, but for the spanned bandwidths and observation frequencies typical of most 
VLBI experiments, this is a good approximation. The assumption is also made that, at a 
given time, the phases of all bins in the band are approximately equal. This implies the 
instrumental channel-to-channel phase offsets must be calibrated with either instrumental 
phase calibration discussed in section 4, and/or a manual phase calibration explained in 
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section 5, for this phase tracking algorithm to work properly (note that manual phase 
calibration will remove linear trends caused by the above media effects; see [5.28], for 
example). 

The first step in the phase tracking algorithm is to form a frequency averaged phase, 
ipj, defined as 

<Pj = argfy^pjj), [13.1] 

t 

where the sum runs over all frequency bins in the band and p t y is the complex data phasor. 
This phasor sum avoids problems near a cycle crossing because it is insensitive to integer 
cycles, and results in a more accurate phase from averaging many frequency bins. The 
problem is now reduced to phase tracking the <pj to obtain an average cycle ambiguity at 
each time point. 

The phase ambiguity of <Pj at time tj is resolved by fitting the previous three points, 
<£>y_ 3) tpj_ 2 and to a line and calculating the expected phase at tj. The integer 

number of cycles needed to move <pj to within cycle of this predicted phase is taken 
to be the cycle ambiguity. This integer ambiguity is then applied to each frequency bin in 
the band being analyzed. Finally, each bin’s cycle at tj is adjusted to be within ±| cycle 
of the average phasor ’s phase, the phase tracked <pj. This process is repeated for each time 
point in the scan resulting in a cycle ambiguity assignment for every integrated phasor. 
Figure 13.2 shows the result of phase tracking for the same channel shown in figure 13.1. 
These phases are now fit with the phase model discussed in the next section. 
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RESIDUAL BIN PHASE. CHANNEL 1 

BSLN 1 BAND 1 XC DATA SCAN 033 AD88 154-3 



Fig. IS. 2 Residual phase as a function of time for the four bins of a single 
channel after resolving cycle ambiguities. The tick marks on the 
vertical axis are full cycles. These tracked phases are next fit 
with least-squares analysis. 
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Section 14 


Least-Squares Fit of Residuals: 
The Phase Model 


Two different models can be used to obtain the final parameter estimates. This section 
describes a fit to the phase model given by [6.2], while the phasor model fit, given by [6.1] 
and [6.2], and identical to the a priori model, is described in the next section. Fitting to 
this model requires the phasor cycle ambiguities to be correctly assigned as described in 
the previous two sections. 

The final parameter estimates are obtained using a linear least-squares fit to the phase 
of the data phasors using the model 


4>\j = <f > o + r T (uJi - w 0 ) + i r Ui{tj - to), 
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for the t'th bin and jth time point. The parameters to be estimated, <f> r Q , T r and f r , are the 
residual phase, BWS delay and phase delay rate, and the superscript r denotes the residual 
values after removing the a priori estimates with the phasor counterrotation. Note that no 
media or instrumental terms are included here; only the three parameters indicated above 
are estimated in this analysis. The effects on these parameters from a linear troposphere 
model are given at the end of this section; the effects from a more general dispersive and 
tropospheric model are calculated in the consider analysis presented at the end of section 
17. Using a least-squares analysis^, the following matrices can be defined: 



is the matrix of parameters, and 


/< M 
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1 . 4 . 3 ] 
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is the matrix of observations (the phase of the integrated phasors) where n implicitly 
includes all time points in all bins. The design matrix, which is the matrix of partial 
derivatives of the model, [14.1], with respect to the parameters, [14.2], is given by 

f 1 Wi — Wo wi(ti — to) ^ 

1 U >2 — U >0 U) 2 (t 2 — to) 


A = 


14.4 


V 1 U n — U ) 0 U>n(t n — to) J 

The weight matrix, W , is the inverse of the observation covariance matrix and is assumed 
to be diagonal for now: 

/l/< 0 ... 0 \ 

0 i/oL ... 0 


w = 
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VO 0 ... 1/alJ 

where aj. is given by equation [2.15c] (off-diagonal elements could be included, for example, 
to model temporal phase correlations due to tropospheric fluctuations). The A T WA matrix 
is given by 


a t wa = 


( 

O j, 

<Pi ~<Pi vf. 

V Y- ( yt - ^o ) 2 v- Ui(ui - U} 0 )(ti - to) 


J 


[14.6] 


7 *i a +i 

. ^ w»(t£— to} y-v w t ( Uj — W 0 )(t,- — to) W t 2 (f, — tp) i 

where the sums run over all time and frequency points in the band. The S-band A r WA, 
for example, can be written in terms of the symbols defined in Appendix C (refer there 
for variables not defined in this text) 

/ So Si — OJqSq 


A t WA = 


Ss — toSi 

Si - u)qS 0 S 2 - 2u 0 Si + OJqSq Sq - t 0 S 2 - ojqS s + w 0 t 0 5i 
V 5 S - t 0 Si S 6 - t 0 S 2 - w 0 S 5 + uj 0 toSi S 9 - 2 t 0 S 6 + t$S 2 J 


14.7] 
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a similar expression is used for other frequency bands. The inverse of this matrix is the 
covariance matrix for the parameters in the X matrix, [14.2]. 

Calculating (A T WA)~ 1 explicitly allows the minimization of the parameter variances 
by appropriate definition of to and ujq : 


( A T WA )~ 1 


1 

Sdo 


S 20 ~ 2u>oSt3 + Wg-StO 
— St3 + U>oStO 

V S 16 — WQ Sti 


— St3 + Wo-Sjo 

St 0 

— Sti 


S 16 — wo^ti \ 
— Sti 

S 02 J 
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where all wo dependence is shown explicitly and all to dependence is contained implicitly 
in the St * variables defined in [C . 10] . Comparing this with the definition of the covariance 


matrix 

{ °lo 

P<f>T^(p 0 

P<t>T a 4>o G T y 



(A t H'A)~ 1 = p$ T o t CT0 o 

ol 

p 7 jG T Oj 

[14.9 


V PtfCfOfa 


of J 



leads to the following: 


Ji _ 

S 20 — 2wo*S't3 + WgS(0 

[14.10a] 

0 

SdO 

2 

S 09 “ 2fo*5o6 + 

[14.106] 

°T = 

Sd 0 

o\ = 

S 02 
Sd 0 

[14.10c] 


— St3 + U>oStO 

[I4.10d] 

P<pr = 

Sd0 o\ 0 ■? 


516 — Wo Sti 

[14.10c] 

II II 

•t- *fc- 

Sdo o\ of 

-Sti 

[14.10/] 

Sdo o\ o\ ' 


These equations lead to the following conclusions: 
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• The phase variance is minimized if to = S$/Si and wo = Sts/Sto- This definition of 
wo simplifies to wo = Si /So for the to definition given. The phase variance is 1/So at 
the minimum. 

• The delay variance is minimized by t 0 = S06/S02 and is equal to S0/S02 at its 
minimum. 

• The rate variance does not depend on to or wo; it is always equal to Soi/Sjo- 

• The phase-delay correlation coefficient is zero for wo = So/Sto ■ This simplifies to 
wo = Si/S 0 for either to = S5/S1 or t 0 = Soo/So2. 

• The phase-rate correlation coefficient is zero for wo = Sw/Sn. This also simplifies 
to wo = Si/ So for to = S5/S1, but wo is undefined for to = Sqq/So2- 

• The delay-rate correlation coefficient is zero for to = Soq/So 2 - 

From these facts, there are only two viable candidates for the definition of t 0 : to = S5/S1 
or t 0 = S06/S02. The latter appears to be the better choice because it both minimizes 
the delay variance and makes the delay-rate correlation vanish. This definition must be 
rejected, however, because this f 0 is often numerically far from the center of the scan; 
sometimes it even appears before the scan start or after the scan’s end. For this reason, 
Fit uses t 0 = S5/S1 as the reference time definition. Because there is only one time 
reference for the scan, to is calculated from the highest frequency band data. Thus, if X 
band is the highest frequency band, 


to = Xs/Xi [14.11J 

(see equation [3.1] and [C.l]). For the band frequency reference. Fit uses 

uox = Xtz/Xto [14-12] 

(see [3.2] and [C.10]), and a similar expression for other bands, because this is the best 
definition in all cases above except for the phase-rate correlation. Even that case simplifies 
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to the optimal value with the t 0 defined in [14.11] for the highest frequency band, and is 
close to optimal for other bands. Fit uses [14.12] to compute w 0 in case the user overrides 
the default to value given by [14.11] or when the band in question is not the highest in 
frequency. For reference, when to and wo are given by [14.11] and [14.12] (as is the case 
for the highest frequency band), {A T WA)~ 1 becomes 


( A 1 W)- 1 = 



( X 0 

0 

0 

1 __ 

0 

Xo(X\X\o - X 5 X l& ) 
X\X d0 

-XoX,6 

XiX d0 


0 

— XnXi6 

X x X d0 

X 02 

X d0 


A 


14.13] 


The matrix (A T WA)~ 1 (A T WI) gives minimum variance, unbiased estimates of the 
parameters in X. A T “WJ is (see [C.2]) 


E 


4>i 





/ 5 0 \ 

a t w = 

— Wo ) 

— 

Sb wo5 a 1 


4 

“ t 0 ) 

2 y 


& 

1 

O 

& 
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Multiplying this by [14.8] gives (see [C.7], [C.8] and [C.10]) 


( S • [5,1 + UJoS s 2 + U>o£()£'a3] \ 

S* • [5 a2 + t 0 S s3 ] , [14.15] 

^ • 5,3 J 


<t>0 -- {St • 5, i) + u> 0 r T 

&d 0 

[14.16a] 

r r = ^-{SfS s2 ) + t 0 f r 

[14.166] 

T r = ~^—{St • Ssz) . 
Zd 0 

[14.16c] 



so that 
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This explicitly shows the delay dependence on the phase delay rate and the phase depen- 
dence on the delay (and thus rate) are as expected. Once the fit values are obtained using 
(A T WA)~ 1 A T W7, they can be added to the a priori values to give the final estimates for 
the residual BWS delay, phase delay rate and phase. Summing these residual values with 
the correlator model values, discussed in section 16, gives the final total observables for 
each band. 

The band amplitude is not determined from the phase fit given in this section, but 
it can be determined separately from a weighted average of all phasor amplitudes in the 
band given by 


A b = N l 




14.17) 


where N L is the number of lags used during the Block II correlation, A\ - is the bandpass- 
corrected phasor amplitude, a\, is the amplitude variance, and i and j run over all 
frequency bins and time points respectively. The primes denote the amplitude and its 
variance after a bandpass correction; A' {j = and o A ^ = o A .B t , where A ij is the 

phasor amplitude, is given by [2.16c] and B { depends only on the bin number within 
each 2 MHz channel, and accounts for one source of instrumental amplitude loss! 9 ). 

After fitting, the data are counterrotated by [14.1] using the estimated residual param- 
eters. The counterrotation is done in the exact manner discussed in section 11 using 


and (f> r {j as the counterrotation phase, instead of <f>A in [11.1]. This final counterrotation is 
done to place the data phasors in their final residual form, making residual plots easy to 

generate by Fit. Figure 14.1 shows the final residual plot corresponding to those shown in 
the previous section. 


♦ 
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RESIDUAL BIN PHASE, CHANNEL 1 

BSLN 1 BAND 1 XC DATA SCAN 033 AD88154-3 



TIME (Sec) 

Fig. 14-1 Residual phase as a function of time for the four bins of a single 
channel after the final fit. The tick marks on the vertical axis are 
full cycles. This final residual plot corresponds to those shown 
in the previous section. 

Consider Analysis: Linear Troposphere Model 

This subsection calculates the effects of a temporally linear troposphere model on the 
parameter estimates; a more thorough treatment of media effects on parameter estimates 
is given in the last subsection of section 17. This simple example illustrates that consider 
errors having the same form as an estimated parameter affect the estimated parameter 
directly. Assume the model phase, 4>\T, is given by 

<t>\™ = 4>o + T(ui - u>o) + - t 0 ) + UiT lt , [14.19] 
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where 


T lt = T l + T* (tj - t 0 ) [14.20] 

is the tropospheric delay, and F* and r* are parameters modeling the constant and linear 
trend of the troposphere. This expression for <£ tJ can be substituted into [14.3] to calculate 

X = {A T W)~ l A T W7, [14.21] 

where X is the matrix of estimated parameters corresponding to X . A shortcut to this 
computation is to notice that the model phase can be rewritten as 

= i<t>o + ft wo) + (r + r*)(w t - — w o) + (t 4- f l )u >{(tj — t 0 )- [14.22] 


Written in this form, it can be seen that [14.3] becomes 7 — Ay where the matrix y is 


( <f > o + T l U } 0 \ 

■ 


14.23] 


Substituting Ay for 7 in [14.21] leads immediately to X = y. In other words, all constant 
terms in the model phase pass through the least-squares formalism into the residual phase 
parameter, all terms (independent of t and j) multiplying w t - - w 0 pass into the delay, and 
all terms multiplying u/ t (* ; -f 0 ) pass into the phase delay rate. Any additional terms in the 
model phase which cannot be put into one of the three categories above (for example, terms 
in the last line of [5.31]) will impact the parameter estimates in a more complicated manner. 
The consider analysis of section 17 will examine such additional terms. In summary, once 
[14.19] is put into the form [14.22], the observable effects can be seen immediately: 

4>o <t>o + o [14.24] 

T — ► T + T* 

• , 

T — *• T + T . 


From this specific example, it can be seen that the troposphere delay and rate enter the 
observables as expected. 
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Section 15 


Least-Squares Fit of Residuals: 

The Phasor Model 

This section describes a fit to the data using the phasor model given by [6.1] and [6.2], 
and is typically used when the SNR is too low to accurately track the phase. The primary 
results of this fit are the parameter variances. Because the model used here is identical to 
that used for the a priori estimate, the estimated parameter means are usually very small, 
usually on the order of 10 femtoseconds for the delay, 10 attoseconds/ second for the phase 
delay rate, and 1 millidegree for the phase (for detected sources); the fit would give zero 
means in the limit that the a priori estimation were done perfectly. The numbers above 
indicate the a priori estimate is usually sufficiently accurate, and improving the search for 
the maximum amplitude phasor, described in section 10, is computationally inefficient. 
The final parameter estimates are obtained with a least-squares fit using the model 

Mi = A b exp i^o + T r (u>i - wo) + i T Ui ( tj — to)) , I 15 - 1 ] 

where Mi is the model phasor and A b , <f> q , t t and f r are to be estimated (see [ 14.1]) . This 
model must be linearized by expanding it about trial estimates which are assumed to be 
close to the final, best estimates. Call the trial parameters A b , T and r and expand the 
model about them to obtain 

Mi ~ Mi + e lk \{A b - A b ) + i A b {<f> T 0 - <£ 0 ) + i(w i - w 0 )(r r - f)A b [15.2] 

H-iw t (*i-fo)(r r -r)A 6 ] 


where 

Mi = A b e'k = A b } (** " *<>)]. [15.3] 

This can be rewritten as 

=2 -v-(Afc — At) +\{t<t> r Q — <J>o) + i(w,- — w 0 )(r r — f) + iw t - (ii -to)(f’ r -?)- [15.4] 

V M t } A b 
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Changing to matrix notation, and using the same formalism as in the previous section, we 
can define the following: 

(A b - Af,\ 

Vo - <£o 

T r -f 

\ t'-'t J 


X = 


Pi 
M, 
P 2 
M 2 

M n 


- 1\ 


( - 1 ) 

- 1 


& 

^L. 

to 

1 

to 

1 

1 — ‘ 


= 

Ab 

-ij 


l Aa. e l(<t>n-<l>n) — lj 

Ab 


A = 


_1_ 

Ab 


f ~T i i(wi-wo) iwi(ti-io)\ 


-j- i i(w 2 — wo) ^2(^2 to) 


V ^ i i(w n - w 0 ) iw n (t B - f 0 ) ) 


W = 


( iMx 


V 0 


0 

i/<4 


0 


0 \ 
0 

V°\J 


!5.5] 
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15.7] 
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where A,- and & are the amplitude and phase of the tth data phasor, p,-, o\, is given by 
equation [2.16c] and i implicitly includes all time-bins points. For complex design and 
observation matrices with real parameters, the covariance matrix is (A^WAr + AjWAj)^ 
and the minimum variance, unbiased estimate of the parameters is 


X = {A^WA r + AJWAi)- 1 (A^WTr + AJWT i ) [15.9] 

where the R and I subscripts refer to the real and imaginary parts. The following matrices 
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are given for reference: 



(AIMITr + AJW?!) 



15.10] 


Because the model is nonlinear, the fit must be done iteratively; the estimates obtained 
for X, which are the differences between the final values and the trial values, must be added 
to the trial values to become the new trial estimates. The process is repeated until the 
elements of X are sufficiently close to zero. Fit iterates this process 10 times; this is more 
than enough for the parameter estimates to converge, while the computational cost of 
performing too many iterations is small. 

Although the amplitude is obtained here as an estimated parameter, it is not passed 
to Modest. The amplitude sent to Modest is given by [14.17] (see reference 9 for more 
details). 

Once the final fit values are obtained, they can be summed with the a priori estimates 
to give the final estimates for the residual BWS delay, phase delay rate and phase. Sum- 
ming these residual values with the correlator model values, discussed in the next section, 


15-3 



gives the final total observables for each band. After fitting, the data are counterrotated 
by the phase model using the values obtained here. The counterrotation is done in the 
same manner given in sections 11 and 14 using equations [11. l] and [14.18]. This final 
counterrotation is done to place the data phasors in their final residual form, making it 
easy for Fit to generate residual plots. 
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Section 16 


Correlator Model Restoration 
to Obtain Total Observables 

The quantity of interest for the fitting package Modest is the retarded baseline delay, 
given by 

’•;(<) = \ •!*(« + <(«))- *(*)]. lie.i] 

and its time derivative, where rj(f) is the location of station i at time t, k is the unit wave 
vector in the direction of propagation, and c is the speed of light; media effects are not 
considered here but are included in the actual correlator models. The retarded baseline 
delay is obtained by first calculating station t’s model delay, r £ = k • r t /c, and its time 
derivative, for each station i. The difference between these quantities for the two stations 
gives the model symmetric baseline delay, and its derivatives. These are combined 

with the a priori and least-squares residuals to give the measured symmetric baseline delay, 
T b , and rate, r b . Finally, corrections for the retarded baseline are applied giving the total 
retarded quantities, t£ and f b , appropriate for Modest. The rest of this section is devoted 
to the details of this calculation. 

The Block II correlator passes to the user the model delays actually applied to the 
station bitstreams at specific reference times. These reference times are at the start, 
finish and at two intermediate times of a user supplied dump interval, where a typical 
interval might be about 20 seconds. The model delays at these four time points (from the 
appropriate dump) are fit with a cubic polynomial to derive the model delay and its first 
three time derivatives at the scan reference time, to. Note that, at present, Fit uses only 
the four points from a single dump interval in the cubic polynomial fit so that to will be 
randomly distributed across the interval; no attempt is made to use points from adjacent 
dumps which may be closer to to than those used (if this were done, to would always lie in 
the center third of the 4-point interval). 

The maximum error introduced in the delay and its derivatives by this method of 
interpolation can be estimated by calculating the exact geometric delay and its derivatives 
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MAXIMUM ERROR FROM MODEL INTERPOLATION 


Quantity 

10 Sec Dump 

20 Sec Dump 

30 Sec Dump 

60 Sec Dump 


4.8 x 10 -18 sec 

5.1 x 10- 17 sec 

2.5 x 10~ 16 sec 

4.0 x 10 -15 sec 


8.9 x 10 -17 

1.1 x 10~ 16 

2.0 x 10~ 16 

1.2 x 10 -15 

7 <7 

6.4 x 10~ 18 sec -1 

2.4 x 10 -17 sec -1 

5.5 x 10 -17 sec -1 

2.2 x 10 -16 sec -1 

F 7 

3.0 x 10 -18 sec -2 

6.0 x 10 -18 sec -2 

9.0 x 10 -18 sec -2 

1.8 x 10~ 17 sec -2 


Table 16.1 The maximum error in the retarded station delay and its deriva- 
tives due to model interpolation for example dump intervals. 

at the scan reference time and differencing these values with the fit values. This calculation 
can be done as a function of baseline orientation to the source position and as a function of 
the scan reference time within the four time points in the dump interval. Table 16.1 shows 
the maximum error found for the delay and its derivatives for several example dump 
intervals. It will be shown below that these model interpolation errors dominate other 
systematic effects in the model restoration procedure. Note that these are the maximum 
interpolation errors expected when calculating the correlator model at a specific reference 
time; they do not account for errors, if any, which are caused by the application of this 
model at the correlator. 

The station model delay (for station »'), r™ r (f), applied to the bitstreams by the Block 
II correlator is a retarded delay designed to cancel significant effects in the fringe phase 
caused by offsetting the bitstreams by the quantized model delay (see equation (7.11) in 
reference 2). It is defined in terms of the unretarded station delay, r£(f), as 

[16.2] 

(see equation (6.2) in reference 2; the equation given here is different, but is the actual 
expression used by the Block II correlator). This retarded station delay should not be 
confused with the retarded baseline correction discussed below. 
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The retarded station delay and its derivatives are given by 


C (1) = <(<) - r”(t)f”(t) - CM 

CM = ?"(<) - C 1 ) r"(() - 3f”(t) f™(t) 

= C‘)+°( 4 )- 


16.3] 


Terms having four or more time derivatives (lumped into 0(4)) have been dropped. To 
calculate the unretarded delay and its derivatives requires that these equations be inverted 
so that r, m (t) and its derivatives are functions of r” r (t) and its derivatives. After some 
algebra, the inverse equations are found to be 


’■.”(<) = CM + CM CM + 2 CM C*M + C’MCM + CM C*M + I 16 - 4 ! 

C*M CM + «C S M CM CM + 
r”(() = CM + CM CM + C 2 M + 2C’M + «CM CM CM + 

C*M CM + 0(4) 


C‘) = CM + CM CM + *CM CM + °( 4 ) 
CM = CM + o(4) 


These expressions, neglecting the fourth-order terms, are used by Fit to obtain the unre- 
tarded station delays and their derivatives. 

To estimate the errors made by using [16.4], equations [16.3] were calculated to fourth 
order in time derivatives and inverted exactly to that order. The contributions from 
the additional fourth-order terms were estimated by substituting the maximum values of 
delay and its derivatives for earth-bound stations, neglecting correlations between these 
maximum values. There are five forth-order unretarded delay terms, the largest of which 
are 48r™ r2 (t) r™ r2 (t) r™ r (t) and 14r™ r (<) r™ r4 (t). The sum of the five terms is always less 
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than about 6.4 x 10 -23 seconds. There are also five additional fourth-order rate terms 
with the largest being 31r™ r (£) i™ r 2 (t) r™ r (f). The sum of the five rate terms is always less 
than about 1.6 x 10 -21 . The additional fourth order terms for r™(<) and r™(f) contribute 
less than 2.6 x 10 -2 ° sec -1 and 2.7 x 10~ 19 sec -2 , respectively. These numbers are all well 
below the interpolation errors shown above in Table 16.1 (and higher order terms should 
be much smaller), indicating that they can be safely ignored. 

The unretarded station model delays are differenced to form the model symmetric 
baseline delay and rate: 


n [16.5] 

where the station indices are shown explicitly. This is combined with the a priori and 
least-squares residual estimates to obtain the measured symmetric baseline delay: 

T b(t) = r?{t) - T a (t) - T r {t), [16.6] 


where r° and r a are the a priori estimates and r r and f r are the least-squares fit residuals. 
The minus signs are necessary to change r° and r r (in the correlator sign convention) into 
the Modest sign convention. Modest requires delays to be in the retarded form shown by 
equation [16.1]. 

Expanding [16.1] and taking derivatives give expressions for the retarded baseline delay 
and its derivatives: 

T b (0 = T b{t) + r r b {t) f S2 (t) + ±r b r2 (t) r, 3 (t) + gr 6 r3 (£) t S2 (t) + 0(4) [16.7] 

Hit) = hit) + T r b {t) f g 2 (t) + r b (t) t S2 {t) + ±T r b 2 {t) r S2 (t) + T b (t) f r b {t) f S2 (t) + 0(4) 

% it) = hit) + T r b {t) T S2 (t) + 2 T r b (t) f S2 {t) + T b (£) t S2 (t) + 0(4) 

Hit) = hit) + 0 ( 4 ), 
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where r a2 (t) is station two’s delay. These expressions can be substituted into themselves 
while dropping fourth-order terms to eliminate the retarded terms on the right. This is 
done iteratively to give 

= 7i(0 +T b (t)T S2 {t)+T b {t)TsS(t) +T b {t)T a2 3 {t) + («)?.,(«) + [16.8] 

\n 2 (t) r S2 (t) f aa (t) + |r 6 3 (0 r S2 (t) + 0(4) 

T r b {t) = h{t) + ffc(0 r, 2 (t) + hit) r 32 2 (t ) + nit) r 32 {t) + T b {t) hit) f S2 {t) + 

2 nit) h 2 {t) r 32 {t) + ^V(t) r S2 {t) + 0(4) 

Hit) = hit) + hit) h 2 it) + nit) r, 2 (i) + 2 Mt) MO + °i 4 ) 

5J(0 = ?i(0+O(4). 

Equations [16.6j can be substituted into the first term right of the equal sign in the delay 
and rate equations of [16.8] to give 

nit) = rrit)-r a (t)-T r (t), [16.9] 

Hit) = T? r it) ~ f a (t) - f r (t), 

where 

M'W = MW + »>(*) MW + n(t) r” 2 W + 'l(f) M S W + jVW MW + [16-10] 

\n\t) M(() MW + gVW MW + o(4) 

^6 mr (0 = Kit) + hit) rZit) + hit) r 3 m 2 2 it) + nit) r™ (t) + r 6 (t) hit) r£ (t) + 

2*1 W M W M W + J*J 2 W M W + 0(4) 

and r"(f), r™(f) and r™(t) model h 2 it), h 2 it) and r,, 2 (i). Equations [16.9] give the retarded 
delay and rate for each band in the final form acceptable to Modest. Note that the a priori 
and residual quantities are applied to the symmetric delay before retardation in equations 


16-5 



[ 16 . 6 ]; this is formally correct to avoid small retardation errors when the measured residuals 
are large, but the difference is usually negligible in practice. 

The total delays and rates obtained here for different frequency bands can be com- 
bined to eliminate dispersive effects caused by the presence of charged particles along the 
observing raypaths. This combination is described in the next section. 
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Section 17 


Charged Particle Calibration 

The radio signal received at the two VLBI stations is corrupted by the presence of 
charged particles along the observed raypaths. These charged particles retard the observed 
single station phase by k a /u>, where k 3 is proportional to the total electron content (TEC) 
along the single raypath and u> is the observation frequency (a conversion factor between 
TEC and k 8 is cr e /27r = 1.3445 x 1CT 7 m 2 cyc 2 /sec, where r e = 2.8179 x 10 -15 m is the 
classical electron radius). The two-station interferometric phase is advanced by 

<t> cv = -k(t)/u, [17.1] 

where k(t) is proportional to the difference in TEC along the two raypaths in the same 
sign convention (either Block II or Modest) as <j> cp . Higher order terms proportional to 
u>~ 2 and u; -3 have been dropped from this equation. These terms can be important at 
lower frequencies, but are usually negligible for observations at S and X-band frequencies 
(2285 MHz and 8440 MHz, respectively) assumed here. The development given below will 
therefore be valid for higher frequency measurements such as X/K-band observations. 

The phase shift shown in [17. l] induces a change in the group delay, phase delay and 
phase delay rate given by 


r cp = 

d<j> cp 

du 

MO 

u/ 2 

[17.2a] 

T cp = 
p ~ 

<^_ = 

M 0 

u> 2 

[17.26] 

f cp = 

P ~ 

II 

n. 

MO 

u, 2 • 

[17.2c] 
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Dual Frequency Calibration 

The single-frequency expressions given by [17.2] can be used to combine the total 
observables from two frequencies into a single observable, free of dispersive effects. The 
S-band group delay in the presence of charged particles is given by 

Ts = f + r cp [17.3] 

= f + k/u 2 Sr , 

where f is the delay free of charged particle effects, k is k(t ) averaged over the scan 
duration, and w| r is an S-band reference frequency for the delay observable, derived below. 
The charged-paxticle corrected delay (or S/X corrected delay), t s / x , can be formed with 
a weighted linear sum of the S and X-band delays. If C ST and C XT are the S and X-band 
delay weights, using equation [17.3] gives: 

T s/x = C*sr T S + C XT t x [ 1 7.4] 

= C Sr (T + k/(jj\ T ) + C XT {f + k/u XT ) 

= r(C ST + C XT ) + k{C st / w\ T + C Xt /u X t ) , 

where t s and t x are the total S and X-band BWS delays. The S/X corrected delay can be 
made equal to f only if the following two equations hold: 


Cst + C X T — 1, 

CsT C X T 

o I O U. 


17 . 5 ] 


( j 


ST 


U 


XT 


This leads to expressions for the delay weights in terms of the reference frequencies: 


C Qt = 


— u; 


ST 


C X T — 


U) 


XT 


w Xr w sr ' 


The other linear combination similar to [17.4] can be used to estimate k: 

ks/X — Cgk T s + C x k T x 

= f (°sk + C- XJk ) + k(c sk /ul T + C xk /u 2 Xr ). 


17.6] 


[17.7] 
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where k s /x is the calculated estimate of k, just as t s /x is the calculated estimate of r. For 
k s /x to be equal to k, the following must hold 


C sk + Cx/t — 0, 

Csk Cxk _ . 
/..2 + ,.,2 - L - 


This leads to 


C S k = 


U 2 ST U 2 XT 


Cxk = 


. ,2 _ , .2 
U/v T Wi' T 


The phase delay rate equivalent of equations [17.4] and [17.7], 

Ts/X — C s j T s + C X T T X 
ks/x = C sk ir s + C xk Tx, 


[17.10] 

[17.111 


and the equivalent of [17.3] lead to similar expressions for the S/X corrected phase delay 
rate and time averaged A;(£) weights: 


Cqr — 


U . 
XT ST 

c _ 

sk ~ 


At o o 

. — (jj L . 

XT ST 

c _ ^If^xf 

X * U 2 .-U 2 .' 

XT ST 


[17.121 


17.13] 


where a and u xf are the S and X-band reference frequencies for the phase delay rate 
observable, derived later. The phase delay equivalents of [17.4], [17.7] and [17.3] lead to 
the phase delay weights, and a second expression for the Jc weights: 


Cstp — 
C s kP = 


U>£ T P 

n „ — 

u XrP 

^Xr p ^StP 

^XtP ~ 

w XrP ~ w srP 


Cxkv — 

iO XT pUg T p 

u xtP ~ U STP 

U XT P ~ W STP 


[17.14] 

[17.15] 


These two expressions are given for completeness, but [17.14] is only needed for phase 
connection (see section 18) and [17.15] is not usually useful. 
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Consider Analysis: 

Dispersive, Tropospheric and Manual Phase Calibration Effects 

The effects on the S/X corrected delay, phase delay rate, and phase observables caused 
by dispersive and tropospheric media effects as well as manual phase calibration processing 
effects are calculated here. The last subsection of section 5 derives the expression for the 
model phase, given by [5.31], which includes these effects. Equation [5.31] must be modified 
to include the counterrotation procedure given in section 11, by subtracting <j > from it. 
The residual model phase just before fitting (section 14) is given by 


4>ij = 4>o~ 4>o + w 0 (r‘ - t u ) - Aw 0 (r c * - 


[17.16] 


+ (r - r* + T* ~ + ^72 “ r°)K - u 0 ) + (f + V - r°)u/,(*j - f 0 ) 


-t 


u> 


_ — k ^ I ic— I U> (t*+t * t* 1 ^4 -u]T tr 

Ui +k Ui + 0 Ui ~ + hi ~ t6i > ~^ + ‘ > ’ 




where Aw 0 is given by [5.32] and u>*, derived later, is [17.29] evaluated for the manual phase 
calibration (MPC) scan. Note that <j>o, t and f represent the total residuals (residuals about 
the correlator model) without media, instrumental or manual phase calibration effects. The 
effect on the total band observables induced by combining many time/frequency points can 
be calculated by substituting the above equation into the linear least-squares formalism, 
given by equations [14.3], [14.14] and [14.15]. It will be shown below that minimizing the 
dispersive effects on the final S/X corrected observables leads to expressions for w| r , 
u/| r p and their X-band counterparts. The phase fit is used in this consider analysis because 
it is used for the high SNR scans; observable variances from low SNR scans are presumably 
dominated by system noise and are insensitive to the small difference in consider errors 
between the phase and phasor models. 

Equation [17.16] is substituted into the A t WT matrix shown in [14.14]. Note that i 
in [14.14] implicitly runs over all time and frequency points while i and j run over these 
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points in [17.16]. The parameter estimates are derived by multiplying this by [14.8] to give 

S ^ • [S sl + U>oS 8 2 + Wo^0^a3] \ 

SI • [ S 32 + toSsz] , [17.17] 

s; ■ s 8 3 

where S £ is similar to S $ in [C.8] except the phase in the expressions for S a , S b and S c in 
[C.2] is given by [17.16]. Expanding gives 

SJ = [*, - <K + wo(f‘ - -r“) - A wo (r‘ - ^ 5 ) - * 8 ] Sv, [17.18| 

+ (r - t* + f‘ -r u + - r'HSv* - w 0 S vl ) 

+ (f + f — T a )(Sv3 ~ toSv 2 ) 

— (fc — k )5v r 5 — ic{Sy 4 — t 0 S Vb ) + ic Sti 
+ CJq S 12 — S t 3 — S t 5 + S tb — S'vfc + Sy T 

and the S terms are given by [C.6] and ]C.12]. Note the last six terms correspond to the 
non-linear part of the dispersive and tropospheric models. Substituting this into [17.17] 
and using [C.9] gives 

4>o + 0o — [<£o — 4>o + — ^ £ *) — Awor*,] — k 5*^+ k*S k *4,— k S^+ S a< j> [17.19] 

r r + r° = [r — r* + (f* - f 4 *)] + £ S kr - k*S k . r + kS kr ~ kS^ r + S aT 

T r + T a = (r + t t ) - k S k T + k*S k *t ~ kS- kj . + ic S ^ + S a j . 

Undefined quantities are given by equations [C.13]. The first group of terms in each 
equation shows the effects on the parameter values from a linear troposphere model, the 
second and third terms give the effects from a constant dispersive model, the fourth and 
fifth terms show the linear dispersive effects and the last term includes the non-linear effects 
from both the troposphere and dispersion models. The signs of the dispersive terms in 
the delay equation of [17.19] are different from those in the other two equations, reflecting 



17-5 



the sign difference in the rightmost terms of equations [17.2], Note the a priori values 
passed through the least-squares analysis correctly, indicating the counterrotation given 
by [11.2] is done correctly (actually, this had to be the case since the design matrix, [14.4], 
is the matrix of partials of [14. l] which has the same form as [11.2]). Also note the delay 
and phase are affected by the difference in the average tropospheric delay between the 
MPC scan and the scan being analyzed, while the phase delay rate is only affected by the 
troposphere rate of the analyzed scan and contains no f U term; this is desired so the time 
derivative of the delay effect is equal to the phase delay rate effect. 

Care is needed in maintaining the two sign conventions coming together here: the 
Block II, station 1 minus station 2, system and the Modest, station 2 minus station 1, 
system. Residual variables, which include those calculated up through section 15 in this 
report, are in the Block II convention. Total variables, those with the correlator model 
restored as described in the last section, are in the Modest convention. All expressions 
in this section are valid in either sign convention as long as all quantities are used in the 
same convention. To avoid considerable switching between systems, all variables in the 
remainder of this section are to be interpreted as being in the Modest sign convention. To 
form total observables here, the left side of [17.19] is added to the correlator model values. 
This is equivalent to [16.6] if the minus signs are absorbed into the residual and a priori 
observables with the troposphere and dispersive modeling, [5.21], by interpreting these to 
be in the Modest convention. Thus, k, ic, k r jt f‘, f, rj r , and r c *„ along with all residual and 
a priori variables, are assumed to be in the Modest convention for the remainder of this 
section. 

Adding the correlator model to equations [17.19] so they can be expressed in terms of 
total S and X-band observables gives: 

4>S = 4>os - <f>o s + wos(f‘ - f**) - Awo s r c * 1(s) - k Sn+k*Sr+ - kS k<f> + kS^ + S a<t> [17.20] 

ts = t- r* (s) + (r* - r ‘ ■*) +k S kT - k*S k . T + kS kr - k\. r + S ar 

t — k S k j- + k S k *j- — kS k j + 1c + S a f 
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4>X = 4>0X~ ^OX"* - W 0x(^ — T* ) — ^^OX^cllX) — £ ^k m <j>~ kXj fc^ + ^ •‘^jfe*^ + X a 0 [17.21) 

T X = T — T* 1(x) + (f* — T**) + k Xk T — k*Xi fc* r + kX kT ~ ic Xj [ , T + Xar 

r x = f + r* — -^fcf+ k* X^j — kXj ci .+ k X^+ X a f, 

where Acuos and Awox are the S and X-band Ac^o defined by equation [5.32], The variable 

on the left and the first term on the right of each equation above are total observables. 
The total observables with a tilde indicate quantities without media, instrumental or MPC 
effects: 


4>os = uo s T™ r + <f>os + 2nN s [17.22] 

4> ox = w o x T b' r + <Aox + 27T N x 
t = r 6 mr + r 

~T = ff r + T, 


(j>o s is, for example, the S-band <j> o, t™ is the correlator model given by [16.10], and N s and 
N x are the integer number of S and X-band cycles needed for phase connection (see next 
section). Equations [17.20] and [17.21] can be substituted into the linear combinations 
such as [17.4] to remove dispersive effects: 


Ts/x - f = (r* - T lt ) - T* - 


w^Ar* 


U) 


XT U ST 


17.23 a] 


+ 


w xr w sr 


{k - k*){ u 2 XT X kT - u} 2 ST S hT ) + k*{u 2 XT ju* 2 T - u 2 ST lu* 2 T ) 


.*2 


.*2i 


+ k{w XT X^ T ttJ ST Sj cT ) ic (u Xt X^ t uj St S + ( oj Xr X aT uj\ r S aT } 


T f 


Ts/X ~T = T + 


W xf " W sr L 


(k k )(uj X j.Xf C T to S j.Sk+) ^{^xtX^ — [17.236] 

+ ic (u> xj .X^ — U s ^S^) + (oJ x j.X a j. — W 5 fSar)j 
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where S k * T = Skr ~ l/ w sr an( * S**r = ■S'jtf ( see [C.13]) are used for simplification, and 


Ar* = r* 


1(X) 


'cl(S) 


[17.24] 


(see [5.30]). The phase delay equivalent of [17.23a] is derived from the phase expressions 
in [17.20] and [17.21] but is not shown here because cycle ambiguities, the N s and N x 
in [17.22], are not known a priori. Typically, the dominant dispersive terms in the phase 
expressions are eliminated directly to resolve cycle ambiguities, making phase delay S/X 
corrections unnecessary. The other linear combination of delays and rates, [17.7] and 
[17.11], can be used to obtain the following: 


t r, ~ <^ \ ( ^ 

*S/X ~ K l 2 ,.2 A, ,*2, .*2 ) + , .2 _ 

Vw xr Xuj xt uj st / w xr 


,2 2 
J XT 
_ .2 

y XT " ST 


A T* + {h-V){S kr -X kr ) [17.25a] 


j _ 

ICs / x ~ ,.2 , ,2 

w xf l 


+ ic(S kT - X kT ) - k {Sy T - x k j + (5 ar - X aT ) 

(k - k*)(S ki - X M ) + k{S k . - X kf ) 

— ic (SjLmj. — Xl,.) — (5 a f — X a j ) 


[17.256] 


In the absence of media and processing effects, the right side of equations [17.23] should 
be zero. Minimizing the largest term in each equation leads to definitions of the reference 
frequencies. The dominant term in [17. 23a] is the one proportional to (k — k*). This term 
can be eliminated by choosing 

w Sr = w xr = l/Xfc r . [17.26] 

Similarly, the k term is the dominant dispersive term in [17.23b] and can be eliminated by 
choosing 

= 1 /Sfr <4 f = 1/X k .. [17.27] 
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Finally, the dominant term in the phase delay equivalent of [17.23a] (not shown) is pro- 
portional to (k — k*) and can be eliminated by setting 

w| T p = 1/ «Sjfc T p w XrP = [17.28] 


These expressions can be written in their more expanded form using [C.13] 


u 


ST 


u 


ST 


2 

U ST P 


Sd o 
t()SdQ 
Sdo 

Sd3 ~ toSd6 

^■’OS’S'dO 

Sd4 — '-“'os Sdb + touosSde 


[17.29] 

[17.30] 

[17.31] 


with similar expressions for X-band. With these definitions and those of [C.14], equations 
[17.23] and [17.25] become 


T S /x — T+[T t - T 1 *) - r* 1(s) - C Xt At* + icR - kT - k R k . r + R aT [17.32a] 

Ts/x = T + t - (k - k*)Rkj- + k R k , , . + R a f [17.326] 

k s /x -k + CsfcAr*, - k*Q k . k + kR- kk - k R^ k + R ak [17.32c] 

k s /x =k + (k- k*)R kj( - k R k . k - R ak [17.32 d] 


These expressions are exact, given their model assumptions. The first two quantities (the 
S/X corrected delay and phase delay rate) are passed to Modest for further analysis. 

The first term on the right side of equations [17.32] shows the uncorrupted observable 
and the terms that follow are consider error terms which can introduce systematic effects 
in the S/X-corrected quantities. The second term on the right of [17. 32a] and [17.32b] 
shows the effects from the average tropospheric delay and rate, and appears to be the 
dominant error term in these expressions at present. An external calibration (for example, 
Water Vapor Radiometer (WVR) data) is needed to eliminate these systematic effects. 
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The last three terms in equations [17.32] can be regarded as data drop-out effects; these 

terms would be zero if were constant in time, for each bin ( a ^ = <7^J. Fit evaluates the 

first two drop-out terms shown (16- byte precision is necessary in this calculation). Because 

they are usually negligibly small, Fit does not use them to correct the S/X observables at 
— — * 

present. The kR kr and k R k , T terms in [17. 32a] are usually on the order of 10 -14 seconds, 
though kR kr may contribute on the order of a picosecond. The (k — k*)R k T term in [17.32b] 

T * 

is typically on the order of 10 -17 , but can be of order 10 -16 ; the Ic R k *j. term is about 

three orders of magnitude smaller. The kR kk and k*Q k *k terms in [17.32c] contribute 

on the order of 10 8 cyc 2 /sec, which typically represents about a 0.1% correction to k s / x \ 

- — 

k Rj c * k is about three orders of magnitude smaller. The ( k — k*)R kk term in [17.32d] is 

T T * 

on the order of 10 2 cyc 2 /sec 2 , which is about a 0.01% correction to k S /x 5 the k R k . k term 
is about an order of magnitude smaller. The last drop-out term in each expression is not 
evaluated by Fit. These terms represent a number of negligible effects, the largest of which 
is probably the one proportional to Syk (see last group of equations in [C. 13]) . Fit does 
not attempt the calculation of Syk from the fcy time series, which is calculated from the 
time series of residual S and X-band delays (Fit calculates for plotting purposes), because 
it is dominated by system noise. 

The remaining error terms in [17.32a], not discussed immediately above, are due to 
differences in the S and X-band clocks. This effect, also discussed in section 5, is written 
as two terms: a constant clock-like term, r c * 1(s) , which is usually absorbed at the Modest 
stage of analysis by a clock parameter, and a term proportional to the S and X-band clock 
difference, Cx- r Ar c *, where C XT is given by [17.6] and A r c * is given by [17.24] and [5.30] 
and assumed to be constant during the experiment. Changes from scan to scan in C X t 
imply a changing delay term which cannot be absorbed by a Modest clock parameter. 
The size of this effect has not been studied thoroughly at this time but may border on 
being significant. Presumably, Ar c * could be estimated in Modest since it is assumed to 
be constant and C XT is known for each scan. 
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The formal errors of the four quantities in [17.32] are given by 


°r S/ x = C sr°r. + C Wr x [17.33d] 

A a/X = c2 sf< + C* xi o} m [17.336] 

*l s/x = C L<+ C >i [17.33c] 

i s/x =^l+C^l [17.33d] 

where a T> and o>, are given by [14. 10b] and [14. 10c], and similar expressions hold for the 
X-band quantities. The formal errors for the delay and phase delay rate, [17. 33a] and 
[17.33b], are also passed to Modest for further analysis. 
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Section 18 


Phase Connection 

This section summarizes the expressions needed for phase connection. Phase connec- 
tion is the process by which the estimated phase, given by [14.16a], [10.2] and [12.4], is 
used to form a phase delay observable. Since only the fractional number of cycles of the 
total phase is measured, the integer part must be determined from other information such 
as the BWS delay and/or model delays. Fit does not attempt phase connection at present 
but passes the necessary information to the Modest binary file where software independent 
of Fit can be used to study cycle ambiguity resolution. 

The residual phase, <f> r 0 , given by [14.16a], and the a priori residual phase, <£“, given 
by [10.2] and [12.4], are determined only to a fraction of a cycle; the true residual phase 
is given by 2 tt N + <f)Q + <f> r Qt where N is an integer. This true residual phase, divided by 
the reference band frequency wo> gives the residual phase delay. If the correlator model 
contains no dispersive effects, it can be applied to the residual phase delay in the same 
manner as the group delay (see [16.9]). The total phase delay for S-band observations, rf, 
can be written 


rg = rr + 


Ms + Ms + 2?r -Ws 

u OS 


= r p + 

— 'Fit(S) ' 


2ttN s 

u>os 


[18.1] 


where Ms and Ms are in the Modest sign convention, t™ is the correlator model delay 
given by [16.10], N s is an unknown integer, and 


T P 

'Fit(S) 


= rr + 


Ms + M 


Os 


WQs 


[18.2] 


is the quantity Fit passes to Modest as the S-band phase delay (S-band is assumed here 
and throughout this section, but similar equations hold for X-band). 
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The S-band phase delay, considering dispersive, tropospheric and manual phase cali- 
bration effects as presented in the last subsection of section 17, is given by (see [17.20]) 


T P - T P - 
1 S — 7 S 


<j> 0 


Os 


+ (T t -T t ')- 


Wos 


Au>0j 

^0S 


^ci(S) kSfcjp "I” ^ &k*TP 


18.3] 


^^ic T p ^ ^k*rP ^arP* 


Here, = </>os/^os is the phase delay without propagation or processing effects, fos is 
given by [17.22] and Aw 0s = coos - ^os ( see [5.32]). The last three terms of [18.3] are 

T T + 

usually negligible and will be dropped and k S^ tP contribute on the order of 1CT 5 

cycles divided by the band reference frequency). The second, fourth and sixth terms of 
this equation depend directly on uq s - If all scans are processed with the same wos, the 
second term will represent a constant, clock-like effect, the fourth term will vanish and the 
sixth term will simplify. Equating [18.3] to [18. l], assuming a constant u 0s for all scans, 
u>0s, and expanding 5* r p and S** T p gives 


p 2nNs _ ~p _ p* ,-t -t*\ {k — k*) 


T Fit (S) "t" 


UOS 


= f p s - r p : + (f £ - f**) - 




stP 


where 


r P* ; ^0 S 
M — _ > 

w Os 


and < jj 2 StP is given by [17.28] and [17.31]. 

Equation [18.4] can be rearranged to give 


[18.4] 


[18.5] 


T P + 
' Fit(S) ' 


{k - k*) 2 ttN s „ 

— ' =T§~ T P * + (f* - f“) 


U) 


+ 


StP 


^Os 


~ r. 


m 


Mod 5 
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where r™ 0<i is the Modest delay estimate which includes clock and tropospheric modeling. 
This can be rewritten as 


r Fit(S) 


+ 


(k - £*) 

U>l T P 


'Mod 


+ 


2wNs 

&0S 


2 ttFs 

<*>0S 


[18.7] 
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where F s is introduced to make a strict equality. If k and k* are approximated by k s /x 
and kg/ x (given by [ 17.32c]), then 


r Fit(S) "b 


fts/X ~ fcg/x) 


U) 


- T m + 
'Mod ' 


StP 


2n N s 

<*>os 


27 tFs 

Oos 


18.8] 


All quantities, except N s , on the left side of this equation are known, so N s can be chosen 
to minimize the right side. This defines F s to be a number of cycles between -1/2 and 
+ 1/2. With N s determined, the right side of [18.6] gives the S-band phase delay. The 
formal error in F s , which can be easily derived from [18.8], [14.10a] and [17.33c], is 



rr 2 2 

°*0 _i ^0 S ( 2 

47T 2 uj T p ' * S/X 
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if 0p s has units of cycles, and has units of radians. A standard check that the assigned 
N s is meaningful is to plot F s for many scans; each F s should be near zero (because the 
Modest clock parameter should remove any constant bias) and their errors should be a 
small fraction of a cycle. 
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Section 19 

Modest Interface 

This section summarizes those quantities passed from Fit to Modest; Table 19.1 sum- 
marizes all such data. As stated at the beginning of this report, the goal of this analysis 
is to derive the charged particle calibrated delay, phase delay rate and phase delay. Fit 
passes three delays and phase delay rates (and their formal errors) to Modest: those for S 
and X band given by [16.9] and the final S/X corrected values given by [17.4] and [17.6] 
for the delay and [17.10] and [17.12] for the phase delay rate. The phase delay for S and 
X band, given by [18.2], is passed but the S/X calibrated phase delay is not because phase 
connection, which Fit does not attempt, is needed first. The amplitude and its formal 
error, given by [14.17], are also passed for both bands. The reference time given by [3.1] 
is added to the scan start time to give the absolute scan reference time. This time, along 
with the reference frequencies given by [3.2], [17.29], [17.30] 
also passed. Finally, the source name, station names, scan 
passed to Modest. 


and [17.31] for both bands, is 
length and scan duration are 
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Quantity 

Reference 

Fortran Storage 

Source Name 


C * 12 

Station Names 


2(C * 8) 

Reference Date 


7*4 

Reference Time 

[3.1]+Scan Start 

R * 8 

Observation Number 


7*2 

Reference Frequencies 

[3.2] [17.29-31] 

8 (72 * 8) 

Scan Length 


R * 4 

Amplitude 

[14.17] [2.16c] 

2(72 * 4) 

Amplitude Error 


2(72 * 4) 

Group Delay 

[16.9] [17.4] 

3(72 * 8) 

Phase Delay Rate 

[16.9] [17.10] 

CO 

to 

*■ 

00 

Phase Delay 

[18.2] 

3(72 * 8) 

Group Delay Error 

[14.10b] 

3(72*4) 

Phase Delay Rate Error 

[14.10c] 

3(72 * 4) 

Phase Delay Error 

[14.10a] 

3(72 * 4) 


Table 19.1 Size and order of data passed from Fit to Modest. The reference 
column lists the relevant equations from the text when applicable. 
An additional 1*2 record descriptor is written to Modest before 
these data. 



Appendix A 

Phase and Amplitude 
Probability Distributions 

The phase and amplitude probability distributions for a complex phasor, along with 
their means and variances, are derived from the Gaussian distributions of the x and y 
phasor components. The phase distribution is derived first by defining the following trans- 
formation and its inverse: 


(x,y) -> (w,wta.n<t>), [Al] 

{4>,w) (arctan2(y, x), x), 


where <f> is the angle to a given x-y point, w is a dummy variable which will be eliminated 
by integration, and the function arctan2(y, x) is defined by 


arctan2(y, x) 


tan 

x > 0 

tan -1 *£ + sign(y) - tt 

x < 0 

sign(y) • | 

x = 0. 


[A.2] 


Figure A.l illustrates the transformation. From this figure, w > 0 implies \<f>\ < | and 
w < 0 implies \<f>\ > the other, blank regions correspond to the transformation <f> = 
arctan2(y, x) + n. 

The (f> probability distribution is obtained by integrating H(<j>,w), the two-dimensional 
probability distribution, over w. The infinitesimals in the two coordinate systems can be 
equated: 


H{4>,w)d<l>dw = f(x)g(y)dxdy, [A3] 

where / and g are the Gaussian x and y distributions. The <f> probability distribution, 
P[d>)i can now be written as 



H(<f>, w) dw 


[A4] 
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Fig. A.l The two-dimensional transformation between the x-y plane and 
the w~4> plane for the function given by [A.l]. The horizon- 
tal lines with arrows show the transformation. The regions in 
the w-<j> plane free of such lines have no corresponding point 
in the x-y plane; these areas correspond to the function <f> = 
atan2(y, x ) + n. 


f | 4 | < | 

1 I-o= H (*.“’) dw M>| 

|jr /<*)»<») i j i*» w<f 

i /foo /(*) »(») W *» W>f 


where | J\ is the Jacobian of the transformation ([A.l]): 
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[A. 5a] 
[A. 56] 


Substituting [A. 5b] and eliminating the absolute value signs in [A. 4] leads to 


P(<f>) = 


sec 2 0 / 0 °° /( w) g(w tan</>) w dw |^| < | 

sec 2 <f> / Q °° /(—to) g(—w tan<£) w dw |d»| > 


[A.6] 


Using the Gaussian distributions for x and y, with means x and y, and variances o\ and 




leads to 


=\2 


— (x — x) 
1 <> _2 

/ *) = 


v/27Ti 


7TCTt 


~(y~y) ; 

, % 1 2cr 2 

9 (y) = at- g y 

V27ra y 


[A.7] 


W) = 


secV r 00 H 

/ exp — - 

27Tt7 z (Tj, 7 0 L 

f 


(±to — x) 2 (itwtan (f> — yY 


2o 2 


+ 


2* 2 


w dw 


2n o x o y cos 2 d> 


—aw 2 ± bw,„ i.„ 
e w aw. 


U.81 


where 


_ 1 / 1 tan 2 <£\ 1 /cos “<p sin"0\ 

2 \cr 2 a 2 ) 2cos 2 </>\ a 2 a 2 / 

_ x ytan</> _ 1 (xcos <j> ys\n<t>\ 


1 /cos 2 (j> sin 2 <^' 


b=^r + 


COS 4 > 


2\a 2 a 2 ) 


[A.9] 


The upper sign in [A. 8] corresponds to \<j>\ < § and the lower sign to \4>\ > f ; this will be the 
case for all such equations below. Making a change of variables from to to r = y/a(w T 
leads to 

Jfe ,2 -c 


P(*) = 


2'Ko x Oy acos 2 <f> J^t 


2 

/ e“ r {r±k')dr 9 


[A. 10] 


where 


k' = 



[All] 
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sign(cos <f>) 


( x cos 0 ysin<ft \ 

V al o?. / 


I cos 2 <£ , sin 2 0 

W _2 H 12 


Define fc = A:'sign(cos </>) = dtzk 1 (using the sign convention noted above) so that 


PW = 




2n a x a y a cos 2 <f> J_/ c 


f°° _ 2 

J e r (r — k) dr 


(for all <f ] . 


[A. 12 


Evaluating the integral and substituting in the expressions for a and c gives the final 
probability distribution for (f> as 




PW = 


k = 




2?r(cr 2 sin 2 </> + cr^ cos 2 <£) 
/ xcos (j> y sin 

1 ' °l °l ' 


1 + y/ir ke k (1 4- erf k )] 


\A 1.13 


y/2 

lcos 2 <f> , sin 2 </> 
[ 2 ' 2 


V a x °y 


Equation [A. 13] was derived for general x and y Gaussian distributions, but as indi- 
cated by equation [2.9], the variances of the x and y components are approximately equal. 
Equation [A. 13] can be simplified using o x — a v = a to give 


. £2 

P{<t>) = [1 + y/n ke* 2 [1 + erf k)], 

£* 7T 


k = 


R 

y/2a 


cos ((f) — (j>). 


[A.14] 


where R and 4> are defined by x = Rcos<j> and y = Rsin<j>. This should be compared to 
equation [2.15a]. 
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The mean of <f> is seen to be unbiased by evaluating 


<*> 



d<f> 



(a + <f>) P(a + 4>) 
a P(a + <j>) da + 


da 



da 


= <t>. 


[A.l 5 


where a = (j) — 4> and the last two integrals are zero and one respectively, because the first 
is odd with even integration limits and the second is the distribution normalization. The 
phase variance is given by 


Var {(j)) = 


[ (<t> ~ P{4>) d<j> 

J — 7T 

= [ % 2 P{<f> + $) d<j> 

— f <t> 2 P{4> + 4>) d4> 

J — 7T 

7T 2 ^ 1 R r l2 1 - R2 °'f<P , Rcos<f>\ 

= — e 2a -I — — — / (f) 2 cos<t>e 2 a ( 1 + erf — 7= — J d<p. 
3 \/27 t a J 0 \ y/ 2 a ) 


[A.16 


Equations [A. 15] and [A.16] are equivalent to [2.15b] and [2.15c]. 


The Amplitude Probability Distribution 


The amplitude probability distribution can be obtained with the transformation 


(x,y) — ► (tw, \/ A 2 - w 2 ) 


[A. 17] 


{A,w) -> (\/x 2 + y 2 ,z), 
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Fig. A. 2 The two-dimensional transformation between the x-y plane and 
the w-A plane for the function given by [A. 17]. The horizontal 
lines with arrows show the transformation. The dashed lines in 
the w-A plane show A — ±u>. 


where A is the amplitude. Figure A. 2 illustrates this transformation. Because both signs 
of y transform to the same A-w point, this transformation is not one-to-one; this is han- 
dled explicitly by a separate integration for each sign of y. The amplitude probability 
distribution can now be written 


/ OO 

H(A,w) dw 

-OO 


-s: 


f( x )g{y) \J\dw + 


f f{x)g{-y)\J\dw, 

J — OO 


[A. 18 ] 


where J, given by [A. 5a] with <f> replaced with A, is A/y/A? — w 2 . Substituting this into 
[A. 18] and simplifying gives 


A fA 

= W* " 


where R is defined above. 


£M 


x w + y v A 2 — w 2 


) 


+ 


(■ 4-191 


exp 


xw — y \/ A 1 — w" 1 


) 


dw 

\fA — tu 2 ’ 
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Using x = R cos <j> and y = R sin <f> as above, and changing variables from iy to (3 with 


w = A cos (3, leads to 




A r / , ££«>•(/?+*)' 


27TC7 2 


— (fl 2 +A 2 ) rir RA 

-e 2 cr 2 / e <^ 2 


1 f ( efcoitf-i) 

Jo 1*' 

— ( J R 2 -j-A 2 ) /*7T — 0 , 

— 1 o L I ( COS Of 

e 2 <J 2 J [e<y 2 + 

f* cos a 

/ e <y da , 

Jo 


+ e & 


[A, 20] 


COS Qt\ , 

e Ida 


where a = /? — <£. Using the integral representation of the Modified Bessel function! 10 ! 


io(x)=±r 

n Jo 


,±x cos/? 


[> 1.21 


leads to the final expression for the amplitude probability distribution (see equation [2.16a]) 


-(tf+A 2 ) 


P(A) = * /o(-p-)- 


[> 1 . 22 ] 


The amplitude mean and variance can be derived in the same manner as the phase 
mean and variance. The amplitude mean is given by 


f 0° 

(A) = AP{A)dA 
J —oo 


[>1.23] 


1 =£■ f°° 

= ~2 e2a \ 

° 2 Jo 


o :z A r ( RA\ f 
A 2 e 2 o- 2 / 0 ^— _j dA. 


The solution to this integral is given in reference 11, giving equation [2.16b]. The amplitude 
variance is calculated using 


/oo 

(A - (A)) 2 P(A) dA 

-OO 

/ OO 

A 2 P(A) dA - {A) 2 

- OO 


= {A 2 )-(A) 2 , 


[A. 24] 



where (A 2 ) is also given in reference 11 as 


(A 2 ) = R 2 + 2a 2 . 


[A. 25] 


Substituting this and [A. 23] into [A. 24] leads to the final expression for the amplitude 
variance, given by [2.16c], 
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Appendix B 


Effects of Time and Frequency FFTs 
on Data Phasors and A Priori Estimates 


This section presents a detailed accounting of the data phasors through the time FFT, 
the rescaling operation and the frequency FFT. The goal is to derive an expression for the 
final, transformed phasors in terms of the original data phasors given by equation [ 10. l] . 

The standard FFT transforms a matrix a(i) into the matrix b(k ), each having m 
phasors. These matrices are related by 


R k ) = Y ®(*) 

i=l 


„2iri 


m 


i,k = 1,2, ••• ,m. 


[B. 1 ] 


By padding the input FFT array with additional zero-amplitude phasors, the output array 
will have more entries, corresponding to a finer discrete sampling of the transform function. 
If m data phasors are padded with zero-amplitude phasors resulting in a total of n phasors, 
the equation relating the data and transformed arrays is 


b(k) = £ a(i) 

i=l 


„2?ri 


(f-i)ffc-i) 

n 


f * = 1, 2, • • • , m 
1 k — 1, 2, • • • , n. 


B. 2] 


From this equation, it can be seen that no approximation or inaccuracy results from 
padding the FFT arrays; the padding simply increases the sampling density of the trans- 
formed phasors. 

The first FFT performed by Fit is over the phasor time series of each bin. If p kl is 
the data phasor for the A:th time point and /th bin, then the transformed array, is 

given by 

B{iJ) = Y PU exp 27ri^ ^ [B.3] 

k=i nt 

where n r is the number of data time points, and nt is the total number of phasors in the 
FFT after padding, and is selected to be the lowest power of two that at least doubles n r . 


B-l 



The time FFT is shown schematically in Figure B.la. The lower numbers in the 
figure indicate the elements of the FFT output array, and the upper numbers represent 
the fringe frequency corresponding to the given array element. The first array element 
corresponds to zero fringe frequency. The second element corresponds to fringe frequency 
A = 1 (cycle) /nt t c Hz, where t c is the number of seconds between points in the time 
series, usually two seconds with the Block II correlator. Each successive array element 
corresponds to a fringe frequency A higher than the previous element. The sampling rate 
of l(cycle)/f c Hz limits the fringe frequency to ±1 (cycle) /2£ c ; actual fringe frequencies 
outside this range will appear in it through aliasing (r in the figure corresponds to this 
highest frequency; r = 1 (cycle)/2£ c ). Thus, array elements in the second half of the array 
should be interpreted as negative fringe frequencies, increasing from the largest negative 
fringe frequency — r Hz to —A Hz. 

Because the fringe frequency in Figure B.la is discontinuous at the midpoint, Fit swaps 
the array halves, as shown in Figure B.lb. Equation [B.3] must be modified to reflect this 
change in array ordering. Replacing i with i — n*/2, the swapped array of transformed 
phasors, F'(i,l), can be written in terms of the original phasors as 

= J2 ? kiexp2iri ~ — ~ — “• ( B - 4 1 

nt 


Fit next interpolates and rescales the transformed phasors. The resulting density is 
half that found in a hypothetical bin having frequency wo- For the purpose at hand, this 
can be broken into two stages: the interpolation, followed by the rescaling. If F"(i,l ) is a 
matrix element after interpolation, then 


F"(i,l) = F\2 «,/) 

n T 

= PM ex P 2n ‘ l 

ifc=l 


[Jc - l)(2t - nt/2 - 

n t 


1 ) 


[B- 5 ] 


This operation is shown schematically in Figure B.lc, where A' = 2 A. The rescaling 
results in the nominal fringe rate being multiplied by lji/uiq. If F'" is the matrix after 
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a 


0 



Fig. B.l Schematic representations of the FFT arrays at various stages 
of the analysis. The lower numbers are the FFT array elements 
and the upper numbers are the corresponding fringe frequency 
(a-c) or BWS delay (d-e). 
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rescaling then 


= £feexp2ri ( *- 1)12, '- n,/2 - 1) ^. 
“ nt w 0 

J k=l 


[B.6] 


Figure B.lc still applies after rescaling if r and A are assumed to be scaled by w//wq. 


Fit next performs the frequency FFT. This FFT not only has zero-amplitude phasors 
padding the end of the array, but also has them between channels when BWS is used to 
increase the bandwidth. If G(i,j ) is a phasor after performing the frequency FFT, then 


G(„ j) = E I) exp 2,1 

U f 


B.7\ 


l=i 

Tlfj 71 r 

= ££?*< exp 2ni 

1=1 k = 1 


r (fc- l)(2t-n t /2-l) w t + (T(l) - l)(j - 1 )] 
L nt wq ny 


where nj, is the number of bins, and nj is the lowest power of two that at least doubles the 
number of possible bins across the observing band. The T(l) array maps the sequential bin 
number to a hypothetical bin number assuming the observing band is packed as densely 
as possible: 

T{1) = Ul ~ Ul + i 5 [B.8] 

Wj 

where w*, is the width of a single bin and wi is the frequency of the first bin. Figure B.ld 
shows the frequency FFT in a similar manner to Figure B.la. Here, 6 = 1 (cycle)/ny wj, 
seconds is the phasor spacing in delay, and ±t is the range of BWS delays, where t = 
1 (cycle)/2wj, seconds. The final operation Fit performs is to swap the array halves in the 
frequency direction to eliminate the discontinuity at the array center. This is the same 
as replacing j with j - n.y/2 in [B.7]. Substituting [B.8] leads to the final form of the 
transformed phasor, T,y, 

fu = y~y'feex P 2Ti[ (t _ + (m-mw-tt/A-ii i 

i,i hi 1 n > “o “i «/ J ’ 

[B.9) 


which is also shown schematically in Figure B.le (compare to [10. 1]). 


B-4 



Calculation of A Priori Estimates 


The maximum amplitude phasor among those obtained from the FFTs corresponds 
to a discrete i-j point near the true maximum phasor. An iterative search in i-j space 
is performed near this maximum FFT phasor to obtain the i-j point giving the true 
maximum amplitude; call this point (xy,xy) and call the phase of the transformed phasor 
at this point <f> m&x = arg T(xj, xj). At the maximum amplitude point, the data phasors 
have the same phase to the best approximation of any point in the delay-rate space. If 
Qijkl is the phase of the exponent in [B.9] and is the phase of the data phasor, pki , 
then on average 

$kl d" ®ZiXjkl = fimax- [£.10] 


Rewrite Q Xi Xjkl 35 


®XiXjkl — 


27 r(fc - l)(2x t - ~ y^-t/2 — 1) qj[ 27r(w/ — u/i)(xy — n//2 — l) 

wo 


_ _ . w/ _ 2 tt, , 

= 27T(/C — 1) Xj H (u>[ — oh )x, 

Wo w 6 v 3 


2n 

tc 


[tk ~ t\) 


— Xi + — (w/ - Wi)xy, 

WO Wfc 3 


£.11] 


where 

2x{ — nt/2 — 1 Xj — n// 2 — 1 

£ Xj = -l , £.12] 

n t 3 rij 1 J 

and [k — 1) = ( t * — t\)jt c are used. To obtain the a priori estimates, replace the data 
phases, $*/, with the model phase given by [6.2]. Substituting [B.ll] and [6.2] into [B.10] 
leads to 


<£max — TU>1 — TUJo + — TCU^o + <f>0 4 ~Xj ( ffc — tj) + — Xj(u>l 

tc U > 0 U)h 3 


[£. 13 ] 


. . 2n 2ir _ 2tt _ 2i r 

= ( r + ~T X i)Ultk + [t - ft 0 — t\Xi H Xj)ui + ( <j> 0 - TU)Q W i X y ) , 

WQtc UJ 0 t c U b J ' u b 33 
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where all k and l dependence has been placed outside the parentheses. For this equation 
to hold for all k and l, the quantities inside the first two sets of parentheses must be zero. 
Thus, the following equations must hold 


. 2tt 

0 = t -| x, 

m ot c 

. 2tt _ 27 r _ 

0 = r — rto t\X{ Xj 

Wo tc Wfc 

, , 2?r _ 

<t > max = <P 0 - TWO W\Xj. 

Ub 


B. 14] 


Solving these equations for r, r and <£o gives 


T = 


27T 


U>ot c 


(to ti)x,' 



= r(to - tl) 



[fl.15] 


. 2tt 

r = -x, 

w otc 


2tt, 27t 

00 Pmw "7 (to tiJXj (w 0 WljXy 

tc Wj 


= <£max + r(w 0 - Wl) + TU>l(t 0 - tl), 


which are also equations [10.2], as desired. 
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Appendix C 

Time and Frequency Summations 
and Expressions 

This section defines time and frequency summations, and combinations of these, which 
are used extensively in the least-squares analysis found in section 14 and the charged par- 
ticle and tropospheric consider analysis given in section 17. Some simplifying expressions 
are also included here which are useful in deriving some of the equations in those sections. 
For the sums that follow, i runs over all bins in the band with bin frequencies w t -, j runs 
over all time points, tj, a t 2 ; is given by [2.15c], and <£ t y is the phase of the data phasors. 
S-band expressions are given here; Xi would be the X-band counterpart to S i, for example. 


(*> i 

= 2 

i Wii 

Srn = S 3 m 

j 

(S) 

IJ l 3 

[C.l 

(5) 1 
S o - S 

So = Y, s 0 
j 

(*> , 

1 j *3 


(s) 

i u 

si = X s ! 

3 

(S) 

_ ST' 

- a ij 
tj l 3 


(5) 2 

si = T^~ 

T 4 

si = y, s i 

i 

(S) 2 

_ V" Uj 

- 4 

t] V 



•S3 = tj Sm 

j 

(S) t- 

_ q 



S 4 = yy tjSg 

(5) t- 



i 

tj a *j 



3 

{s) u t- 

- 4 

ty l 3 



sq = tj-s^ 

j 

(5) ^ 
w * q 

a 2 
tj *J 
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«7 = E‘i S 0 

3 

Se = E‘i S ! 

3 

S» = E‘ ? S 2 


(S) ±2 

= Y- L 

^4 

tj l 3 

< s ) ,,f 2 

= E^ 


(5) 


= E 






2,2 

w » t i 




The following expressions contain the interferometer phase: 


(5) A 

s.^Esi 

3 

{s) A 
_ -ST- <Pij 

~ ^ 4 

tj 

(S) A 

qj - V" 'Vij 

4 = V 4 

* = E s / 

3 

(*) , 
v - ^ w t<Pty 

4 

tj 


& 

111 

-M 

Vwi. 

Co 

(S) , A 

~ ^ a?- 

»i °'J 


[C.2] 


The following combinations are defined in terms of the summations given by [C.l]: 

S 02 = ^0^2 - S\ [C.3] 

s 06 = SqS§ — S\Sc, 

5ie = S\Sq — S2S5 

S 09 = SoSg - S| 

Siq = 5i5g — S 5 S 6 

$29 = -S'259 — S|. 
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More complicated combinations are given by 




So 

Si 

s 5 

Sdo 

= +Det 

S 1 

s 2 

S 6 



s 5 

S 6 

S9 



S3 

S 4 

s 7 

S<n 

Q 

+ 

III 

Si 

S 2 

S6 



s 5 

s 6 

S 9 



So 

Si 

S 5 

Sd 2 

= — Det 

S 3 

S4 

S 7 



S5 

s 6 

S 9 



So 

Si 

s 5 

Sdz 

= "hDet 

Si 

S 2 

S 6 



S3 

S 4 

S 7 



Sm 

So 

S4 

$d 4 

0 > 

Q 

+ 

III 

Si 

S 2 

s 6 



S5 

s 6 

S9 



So 

Si 

S 5 

SdS 

= —Det 

Sm 

So 

S 4 



s$ 

S 6 

S 9 



So 

Si 

s 5 

$d 6 

= +Det 

Si 

S 2 

S 6 



Sm 

So 

S4 


= S5S16 — SqSoq + S9S02 


= S3S29 — S4S19 + 5 7 5 i 6 


= S3S19 — S4S09 + SjSqq 


= S 3 Si6 ~ S4S06 + 5 7 5 02 


= S m S29 — S0S19 + S4S 16 


= S m S 19 — SqSq9 + S4S0Q 


= S m SiQ — SqSqq + S4S02 


CA 


When each of the quantities above is normalized to the same units by using the appropriate 
factors of ujq and to, it is found that 


0 SdO Sd.1 

2 — « 7 

Wq tQ^O 



SdA 
U 0 


~ 25^5 


[C. 5 ] 


and S d6 , in these normalized units, is five to six orders of magnitude smaller than the 
above variables. With these definitions, all the above quantities are positive except S d 6 , 
which can have either sign. 
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The following vectors can be defined 


5yi = {Sq,Si,Ss) 

Sv 2 = {S U S 2 ,S 6 ) 

Sy 3 = {S5,Sq, 5g) 

Sv 4 = (<53, 5 4 , 57) 

Sv 5 — (5 m , 5o,54) 

551 = (529, — 5i9, 5i6) 
Ss2 = (—5x9, 5 09 , — 5 o6) 
553 = (5i6,-5 0 6,5 02 ) 

S<f> = (5 a , 5j, 5 C ). 

With these definitions, the following hold: 


5(io 

0 
0 

■S'di 

Definitions containing to are: 


Ss2 • 5yi = 0 
Ss 2 ' Sv 2 = 5^0 
5 S2 • 5y 3 — 0 
Ss 2 • 5v4 - —Sd 2 

5 S2 • 5^5 = —5^5 

5xi = 5 o6 — to 5 02 
5x2 = 5o9 — to 5oe 
5x3 = 5x9 — to 5x6 
5xo — 5x2 — to 5xx 

= 5q9 — 2 Ly, S t jQ -f fg 5 02 . 


[C.6] 


[C.7] 

[c.8] 

553 • 5yx = 0 [C7.9] 

Ss 3 • Sv 2 — 0 

5^3 • 5y3 = 5do 
553 • 5^4 = Sd3 

Ss 3 • 5^5 = Sdo- 

[C. 10 ] 


551 • 5yi = 
Ss 1 • Sv2 - 
Ss 1 • Sv 3 = 
Ssi • Sv 4 = 
Ss 1 • Sv 5 = 
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The following identities are useful for simplifying expressions: 


The following 
definitions): 


SqiSqq — Sqq = SoSjo 


[C.ll 


S02S19 — «S’o 6 ‘S'i 6 = SiSdQ 
S 02 S 29 — = 

SoeSig — SoqSig = SsSdo 
‘S' 06^29 — •S’leSig = SeSdo 


S09S29 ~ S- 


19 


— Sg Sjq . 


are defined for the consider analysis in section 17 (see [5.27] for the t* 


Sn 


St 2 
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Ui°ij 


( 5 ) .* ( s ) t t * \ 
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E hj V" _ 2 *_ V" t J t2> | 
, , -2 ’ _2 ’ _2 
ij a ij ' 
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,($) , (5) 2 (5) 2 
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' ij °'J ij °'i ij a 'i 7 


(5) 
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<SVt = 


(5) 


(5) 


( (5) ur 

E^.E^-.E^) 

. . cj t cr,M <7,- / 

ij x V 17 V •* 17 7 


cj t -cr, v 

4 4 * tj 

(S) OJ T tf 

0 } 

13 








( 5 ) (jJ?T tr ^ , 2 t T tr \ 

E u t r 3 l 3 r i \ 

^ v tj 7 
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Several simplifying expressions are defined as: 


S k <fi = {Sd4 ~ ^O'S'ds + UotoS^/Sdo 


C. 13] 


Skr = {Sd5 ~ toSd6)/SdO 


Skf = S dQ /Sdo 


S krP = S kl p/ojo 


s k *<p = S kt f, + Awo / ^sr 

Sk*T = S kT - 1/ W ST 
S k *T = S k r 
Sk‘TP ~ S k ^/U)Q 


S k <p = 1(^1 — w o5(f2 + OJ 0 t 0 S d s) - to(Sd4 - VoSfa + Wofo<S'i6)]/‘S' (i o 
S kr = [(^2 — toSd3) — to(Sds — toSj 6 )]/SdO 
S k r — (Sd3 — to Sde)/ s d0 

S krP “ S k<t >/ U 0 

S'k*<p S ‘ {&S1 + Wo*?52 + WotoSgz) / S d o 
Sjc * T = —Sti • (&S2 + ^oSsaJl/^o 
^i*f = (‘Sfl ' ^53)/^d0 
^k*T p = S k y/ W ° 

^ = ( w 0^2 - ^3 - 5 * 5 + 5 < 6 - Sy* + 5 yy) • (55! + U} 0 Ss 2 + OJot 0 Sss) / S d o 

S aT = (<v$St 2 — Sts — Sts + Ste — Sy k + 5 vv) • ( Ss 2 + toSss)/S ( io 

S a f = (u>QSt2 — Sts — Sts + StQ — Sv k + Sy T ) ■ 553/5^0 

Sar p = S a< f,/ujQ. 
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The following are simplifying ratios: 


R kr = 


R; m = 

k*T 

Rar = 


x kT s kr - S kT X 


hr 


Xkr ~ 

XfcrS^T ~ S kT X' k * r 
Xkr ~~ &kT 

XkrSar ~ 

XjiT ~ $kT 


Rkr = 


Rl* ♦ = 

k* t 


Rar = 


X k j.S k r Sj'j.Xkj- 


kjH 

x kf - s kf 


Xl.Sl .. - S;.X L .. 

kr k*T kr k*r 

X kr - Sic r 

x kf - Skr 


R . krP — 


R'k'rP ~ 


RarP = 


X H S^ - s kt x : 


k(fi 


X k( p - S k( p 
x k<pS k * ( p — S k ,x-^ 
Xh* Rkip 

Xk<pSa<t> ^k<pX a< if\ 

X k <t> - Sk<fi ) 


^JfcJfc — 


~?kr\ 

SkrJ 


.(*kZ 

\x t ,- 

R . = ( ^Jc^T ~ ^ k *T 

k ' k ~\X kT - S kT 

Rak = ^ 


X (IT Sn 


,x kT - S kT 


R kk — 


Rl. 


k*k 


R l = 

ak 


_ ( x k r-s kf \ 

\ X kf-S k J 

( X±LZlkt \ 

\ x k,-s k J 

( Xgi-SgA 

V x k f — s k f J 


R kkP ~ 


R-k*kP 


RakP = 


= ( 

\X k< p - S k J 

( X k*p ~ S k'4> 
\ X k< f, S kl p 

( Xatp Ra<j> \ 

X k p - SkJ, J 


C. 14 ] 
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Appendix D 

Variable Definitions 


The variables used in this report, except for those found in the appendices, are given 
below along with their definitions. The equation number next to each variable gives the 
definition or first occurrence of the given variable; a — or + preceding the equation number 
indicates the variable is defined just before or after the indicated equation, respectively. 
Variables are ordered alphabetically by the primary symbol, with Roman characters first, 
followed by Greek symbols. 

a, —[10.4] Coefficient i (t between 1 and 5) for 2-D parabolic fit 
A —[2.15] Generic phasor amplitude 
Ai [15.6] Phasor amplitude of time-bin point i 
A{j +[14.17] Phasor amplitude of bin i at time point j 
Ajj [14.17] Bandpass corrected phasor amplitude of bin t at time point j 
(A) [2.16b] Mean phasor amplitude 

A^ [2-19] Station i’s tone phasor amplitude 
A [14.4] Design matrix for least-squares analysis 

Aft [6.1] Band amplitude 

Aft [15-2] Trial band amplitude 

Af [7.3] Model FFT phasor amplitude 

A m [10.4] FFT phasor amplitude 

A t [7.1] Test amplitude 


A t 

[7.2] 

Test amplitude for phase rate k 

Aj 

[15.9] 

Imaginary part of design matrix A 

Ar 

[15.9] 

Real part of design matrix A 
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b 


[2.1] Bin index 
B [2.8] Average blanking counts: B = ( B x + B y )/2 

B —[2.12] Average Block II blanking counts: B = ( B x + B y )/2 

B{ -{-[14. 17] Bandpass amplitude correction for channel bin i 
B x [2.5] Real (cosine) blanking counts 

B y [2.5] Imaginary (sine) blanking counts 

B x —[2.12] Real (cosine) Block II blanking counts 

B y —[2.12] Imaginary (sine) Block II blanking counts 

Bj [12.2] Average integrated Block II blanking counts 

c [16. l] Speed of light = 2.9979 x 10 8 m/sec 

C [2.13] Like-sign cosine (real) Block II counts in bin mode 

Ci [12.1] Integrated like-sign cosine (real) Block II counts in bin mode 

Ct { [4.1] Like-sign tone cosine (real) Block II counts for station i 
Csr [17.4] S-band delay weight for dual frequency calibration 

C X t [17-4] X-band delay weight for dual frequency calibration 

Csf [17.10] S-band phase delay rate weight for dual frequency calibration 

Cxt [17.10] X-band phase delay rate weight for dual frequency calibration 

C ST p [17.14] S-band phase delay weight for dual frequency calibration 

C X tp [17.14] X-band phase delay weight for dual frequency calibration 

Csk [17-7] S-band Total Electron Content (TEC) weight 

C X k [17*7] X-band TEC weight 

C s j. [17. 11] S-band TEC rate weight 

[17.11] X-band TEC rate weight 
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C S kP [17.15] S-band TEC weight (from phase delay) 

C X kP [17.15] X-band TEC weight (from phase delay) 

D{ [5.4] Weighted phase average for bin i 

Di [5.7] Weighted phase average for bin i after model subtraction 

E( [5.4] Time weighted phase average for bin i 

E{ [5.7] Time weighted phase average for bin i after model subtraction 

7 [14.3] Matrix of observations for least-squares analysis 

Tr [15.9] Real part of observation matrix 7 

7i [15.9] Imaginary part of observation matrix 7 

Fs [18.7] Fractional number of cycles to connect phase 

F{ [5-24] Special summation involving residual TEC 

[5.24] Special summation involving residual TEC 
h [8.2] Ratio used for phasor interpolation 

* Index, often used to index bin numbers 

/ [2-11] Invalid counts 

I —[2.12] Block II invalid counts 

1 1 [12.1] Block II integrated invalid counts 

hi [2-18] Block II tone invalid counts for station i 

Iq [2.16] Modified Bessel function of order zero 

Ii [2-16] Modified Bessel function of order one 

3 Index, often used to index time points 

kj [5.20] Proportional to baseline TEC at time tj 

k [5.21] Average kj-. proportional to scan TEC 
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ic [5.21] Averaged kj rate: proportional to scan TEC rate 

k T j [5.21] Residual of kj after subtracting linear model 

k* [5.26] Average kj for the MPC scan 

k [5.26] Averaged kj rate for the MPC scan 

k[t ) [17. l] Proportional to baseline TEC at time t 

k(t) [17. 2c] Proportional to baseline TEC rate at time t 

k s — [l7.l] Proportional to station Total Electron Content (TEC) 
k s /x [17.7] Estimate of k using S and X-band BWS delays 

ks/x [17.11] Estimate of k using S and X-band phase delay rates 

Jc [16.1] Unit vector in direction of wave propagation 

l [2.1] Lag index 

L x [2.5] Real (cosine) like-sign counts 

L y [2.5] Imaginary (sine) like-sign counts 

L x [2.12] Real (cosine) Block II like-sign counts 

L y [2.12] Imaginary (sine) Block II like-sign counts 

M +[2.12] Block II normalization factor 

Mi [5.24] Special summation involving residual tropospheric delay 
M* t . -+[4.1] Block II tone normalization factor for station i 

Mi [15.1] Model phasor Ifor tth time-bin point 

Mi [15.3] Trial model phasor for tth time-bin point 

rib [10.1] Number of bins in band 

rif [10.1] Number of BWS delay points in FFT 

n r [10. l] Number of time points 
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nf [ 10. l] Number of phase delay rate points in FFT 
N —[2.11] Number of bits summed per dump 

N{ [5.24] Special summation involving residual tropospheric delay 
Nl [2.1] Number of lags per channel 

Nj [12.1] Number of bits summed per integrated dump 

N s [17.22] Number of S-band cycles needed for phase connection 

Nx [17.22] Number of X-band cycles needed for phase connection 

O [16.3] Combination of terms of this order 

P{A) [2.16a] Phasor amplitude probability distribution 

P{<j>) [2.15a] Phasor phase probability distribution 

P ™ [6.1] Phasor model for bin i at time point j 

r e —[17.1] Classical electron radius = 2.8179 x 10 -15 m 

r,(t) [16.1] Location of station i at time t 

R -[2.16] Location of phasor probability distribution’s center in polar coordinates: 

(R,4>«) 

Si [5.4] Special MPC sum for bin i 

S [2.13] Like-sign sine (imaginary) Block II counts in bin mode 

5/ [12.1] Integrated like-sign sine (imaginary) Block II counts in bin mode 

S ti [4.1] Like-sign tone sine (imaginary) Block II counts for station t 
t Time 

to [3.1] Scan reference time 

<i [10.2] Time tag of first time point 

tj [3.1] Time tag of jth time point 

tj [2.3] Time tag of j'th bit 


D-5 



tc 

[10.2] 

t'o 

[12.3] 

C 

[5.27] 

Ti 

[5.4] 

Tij 

[10.1] 

u x 

[2.5] 

Uy 

[2.5] 

Ui 

[5.4] 

V 

[2.3] 

Vi 

[12.2] 

Vt, 

[2.17] 

w 

[14.5] 

X i 

[10.3] 

Xj 

[10.3] 

Xi 

[10.3] 

Xj 

[10.3] 

X 

[14.2] 

X 

[14.2] 

z x 

+ [2.12] 

Zi 

[5.11] 

Mt b j) 

[2.3] 

6 

[2.15c] 


Time between adjacent time points 
Recalculated scan reference time 

Special time expression number n for bin i used in MPC analysis 
Special MPC sum for bin i 

Transformed (by 2-D FFT) phasor for fth phase delay rate point and j'th 
BWS delay point 

Real (cosine) unlike-sign counts 
Imaginary (sine) unlike-sign counts 
Special MPC sum for bin t 
Valid counts 
Integrated valid counts 
Valid tone counts for station i 

Weight matrix for least-squares analysis (inverse of observation covariance) 

Real variable corresponding to integer phasor index i 

Real variable corresponding to integer phasor index j 

Linear function of x t 

Linear function of Xj 

Parameter matrix for least-squares analysis 

Matrix of minimum variance, unbiased estimates of parameters in X for 
least-squares analysis 

DC bias for Block II cosine accumulator 

Combination of special MPC sums for bin i 

Lag l's single-bit correlation coefficient at bit time t b - 

Dummy phase variable used for integration 
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V 

[4.3] 

Pb{t) 

[2.1] 

Pl{ t) 

[2.1] 

pu (0 

[2.17] 

p 

[8.2] 

Px 

[2.13] 

Py 

[2.13] 

P 1 

[8.2] 

P 2 

[8.2] 

PI 

[12.1] 

Pi] 

[10.1] 

Pti 

[11.1] 

P CC 

[14.18] 

P<pr 

[14-I0d] 

P<t>T 

[l4.10e] 

Prf 

[14. lOf ] 

a 2 

[2.9] 


[12.2] 


[2.17] 

a l 

[2.15c] 

a \i 

[14.5] 


[3.1] 


Calibrated variance scale factor 

Frequency-domain complex correlation coefficient for bin b at time t 
Lag-domain complex correlation coefficient for lag l at time t 
Tone complex correlation coefficient for station i at time t 
Generic data phasor 

X-component (also real-component or cosine-component) of pj 
Y-component (also imaginary-component or sine-component) of pj, 
Generic data phasor 
Generic data phasor 

Integrated complex correlation coefficient 

Complex correlation coefficient pj for bin t at time point j 

Counterrotated complex correlation coefficient for bin i at time point j 
(after a priori model is removed) 

Counterrotated complex correlation coefficient for bin i at time point j 
(after residual model is removed) 

Phase-BWS delay correlation coefficient 

Phase-phase delay rate correlation coefficient 

BWS delay-BWS delay correlation coefficient 

Variance of either component of complex correlation coefficient 

Variance of either integrated phasor component 

Variance of either tone phasor component 

Variance of phasor phase 

Variance of phasor phase for tth time-frequency point 
Variance of phasor phase for bin i at time point j 
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JrtkO ?rno 


[2.16c] 

a ly +[ 14 - 7 1 

< [14-7] 

<4 [ 14 ' 17 1 

ct 5 0 [14.10a] 
a? [14.10b] 
a? [14. 10c] 

a h,x l 17 ' 33 “l 

< lx l 17 - 33b l 

E 17.33c] 

S/X 

[17.33d] 

S/X 

°l s [18.9] 
t [5.13] 
T a [6.2] 
T al [12.3] 
T f [14.1] 
r b [16.6] 
T' h [16.1] 
r"* [16.6] 

r 6 mr [16.10] 
t.. [16.7] 

C (16-2] 

C I 16 - 2 ! 


Variance of phasor amplitude 

Variance of phcisor amplitude for bin i at time point j 

Variance of bandpass corrected phasor amplitude for bin * at time point j 

Variance of band amplitude 

Variance of <f > o parameter 

Variance of r parameter 

Variance of r parameter 

Variance of S/X BWS delay 

Variance of S/X phase delay rate 

Variance of S/X corrected k 

Variance of S/X corrected k 
Variance of F s 
BWS delay for generic band 
A priori BWS delay 

A priori BWS delay corresponding to recalculated to and wo 

Residual BWS delay 

Symmetric baseline delay 

Retarded baseline delay 

Model symmetric baseline delay 

Model retarded baseline delay 

Symmetric station delay for station i 

Model symmetric station delay for station i 

Model retarded station delay for station i 
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t s [17.3] S-band BWS delay 

r x [17.4] X-band BWS delay 

f [17.3] BWS delay without media, instrumental or MPC effects 

f* [5.21] Scan averaged tropospheric delay 

f u [5.26] Scan averaged tropospheric delay for MPC scan 

f [15.2] Trial BWS delay 

t [17.22] Total BWS delay without media, instrumental or MPC effects 
t* [5.19] BWS delay from MPC scan 

r cp [17. 2a] BWS delay change induced by charged particle effects 

Tp P [17.2b] Phase delay change induced by charged particle effects 
[18. l] S-band phase delay 

[18.3] S-band phase delay without media, instrumental, or MPC effects 
T s/x [17.4] S/X corrected BWS delay 

rf [5.20] Tropospheric delay at time point j 

T j r [5-21] Residual tropospheric delay at time point j after subtracting linear model 

T lt [14.20] Linear tropospheric model delay 

r BSA [5.12] Band-averaged bitstream alignment (BSA) delay 

r BSA [5.26] Band-averaged bitstream alignment (BSA) delay for the MPC scan 

f BSA [5-28] Band-averaged bitstream alignment (BSA) delay for the MPC scan without 
dispersive effects 

r c*i [5.30] Clock error for given band, assumed constant for all scans including MPC 
scan 

r c*(s> [17.20] S-band clock error, assumed constant for all scans including MPC scan 

T c*(X) [17.21] X-band clock error, assumed constant for all scans including MPC scan 
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T P * 

T c\{S) 

[18.5] 

S-band clock-like term used in phase connection, assumed constant for all 



scans including MPC scan 

T* 

cor 

[5.28] 

Small delay correction term 

f m 

'Mod 

[18.6] 

Modest model delay 

fP 

' Fit( S) 

[18.2] 

S-band phase delay calculated by Fit 

Ar* 

[17.24] 

Difference between S and X-band clock errors 

j- «ipc 

[5.1] 

Band phase delay rate estimate for MPC analysis 

T 

[5.13] 

Phase delay rate 

T a 

[6.2] 

A priori phase delay rate 

T a> 

[12.3] 

A priori phase delay rate corresponding to recalculated to and wo 

T r 

[14.1] 

Residual phase delay rate 

T S 

[17.10] 

S-band phase delay rate 

Tx 

[17.10] 

X-band phase delay rate 

Ts/X 

[17.10] 

S/X corrected phase delay rate 

ir 

r 

[15.2] 

Trial phase delay rate 

T 

[17.22] 

Total phase delay rate without media, instrumental or MPC effects 

T 

[17.22] 

Phase delay rate without media, instrumental or MPC effects 

-t 

T 

[5.21] 

Scan averaged tropospheric delay rate 

T t* 

T 

+ [17.19] 

Scan averaged tropospheric delay rate for the MPC scan 

"p p 

[17.2c] 

Phase delay rate change induced by charged particle effects 

4> 

-[2.15] 

Generic phasor phase 

(t) 

[2.15b] 

Mean phasor phase 


[4.2] 

Tone phasor phase for station t 
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0* 

-[ 2 . 16 ] 

fiij 

[ 5 . 2 ] 

4>ij 

[ 5 . 6 ] 


[ 5 . 18 ] 

V, 

[ 6 . 2 ] 

<p a .' 

v '*7 

[ 12 . 3 ] 

6 r .. 

[ 14 . 1 ] 

00 

[ 5 . 13 ] 

Vo 

[ 6 . 2 ] 

Vo' 

[ 12 . 3 ] 

Vo 

[ 14 . 1 ] 

Vo 

[ 5 . 19 ] 

00 

[ 15 . 2 ] 

0t 

[ 14 . 3 ] 

*i 

[ 5 . 12 ] 


[ 5 . 6 ] 

V 

[ 15 . 2 ] 

00 S 

[ 17 . 22 ] 

0OX 

[ 17 . 22 ] 

Vos 

[ 17 . 20 ] 

Vox 

[ 17 . 21 ] 

00 S 

[ 17 . 20 ] 


Location of phasor probability distribution’s center in polar coordinates: 

[R,V) 

Phase for bin i at time point j 

Counterrotated phase for MPC calibration for bin i at time point j 
MPC corrected phase for bin i at time point j 
A priori phase model for bin i at time point j 

A priori phase model, corresponding to recalculated to and wo, for bin i at 
time point j 

Residual phase model for bin i at time point j 
Band phase 

» 

A priori band phase 

A priori band phase corresponding to recalculated to and wo 
Residual band phase 
Band phase for MPC scan 
Trial band phase 

Phasor phase for ith time-frequency point 

Bin phase offsets for bin i; output of MPC analysis 

Model residual fringe phase for bin i 

Trial phasor phase for tth time-frequency point 

Total S-band phase 

Total X-band phase 

Total S-band phase from MPC scan 

Total X-band phase from MPC scan 

Total S-band phase without media, instrumental or MPC effects 
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<fox [17.21] Total X-band phase without media, instrumental or MPC effects 

[5.20] Model phase for bin i at time point j assuming instrumental, dispersive 
and tropospheric phase offsets 

(f > +" [5.13] Model phase for bin * at time point j assuming only instrumental phase 

offsets 

<f> ,m — [7.1] Linear model phase 

[14.19] Model phase for bin t at time point j 
<p max [10.2] Phase of the maximum amplitude phasor after FFTs 
[5.1] MPC model phase for bin i at time point j 
^ [5.1] MPC bin i reference phase 

(j> cp [17.1] Phase change induced by charged particle effects 
[5.13] Instrumental phase offset for bin i 
<f>! m [2.3] Model fringe phase corresponding to bit time ty 

(j> [8.1] Band phase rate = uqt 

<Pj [13. l] Bin-averaged phase for phase tracking at time j 

X 2 [5.2] Chi-squared for MPC analysis 

4>i [8.1] Rescaling scale factor 

( j — [17. l] Generic observation frequency 

u>i [3.1] Observation frequency of bin i 

u)q +[3.2] Band reference frequency 

[ 10. l] Observation frequency of lowest frequency bin in band 
cjj, [10.2] Bandwidth of single bin (nominally 0.5 Mhz) 
u>q [12.3] Recalculated band reference frequency 

u^os +[3-2] Band reference frequency for S-band 
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ojqx [3.2] Band reference frequency for X-band 

Wq [5.26] Band reference frequency for MPC scan 

Awo [5.32] Difference between scan and MPC scan reference frequency = u>o — Wg 

Au>o s [17.20] S-Band Au>o 

(Do 5 —[18.4] Common S-Band reference frequency for all scans 

wj [5.6] Model residual fringe frequency for bin i 
uj — - [8.1] Residual fringe frequency for bin i 

uj* [5.28] MPC scan’s u> ST (without the band subscript) for a generic band 
w sr [17.3] S-band reference frequency for S/X BWS delay calculation 

w X T [17-4] X-band reference frequency for S/X BWS delay calculation 

w S f [17.12] S-band reference frequency for S/X phase delay rate calculation 

u; X f [17.12] X-band reference frequency for S/X phase delay rate calculation 

uj St i> [17-14] S-band reference frequency for S/X phase delay calculation 

u> Xt p [17.14] X-band reference frequency for S/X phase delay calculation 

fi — [7.1] Trial fringe rate 

fit [7.2] Trial fringe rate number k 
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